# -----------------------------------------------------------------------
#    Title: Too Tired to Vote: A Multi-National Comparison of Election Turnout with Sleep Preferences and Behaviors
#    Authors: 
#    Date: 03 31 2022 
# -----------------------------------------------------------------------

#### Cleaning the Working Environment ####
rm(list = ls(all=T))
gc()

#### Getting the Working Directory ####
getwd()

#### Creating Empty Folder for Tables in the Working Directory ####
dir.create("Tables")

#### Loading Packages ####
library(haven)
library(sjPlot)
library(nnet)
library(stargazer)
library(tidyverse)
library(viridisLite)
library(viridis)
library(forcats)
library(data.table)
library(openxlsx)
library(psych)
library(lme4)
library(lattice)
library(ggeffects)
library(sjlabelled)
library(naniar)
library(dotwhisker)
library(broom)
library(RColorBrewer)
#library(merTools)
library(ggplot2)
library(openxlsx)
library(apaTables)
library(readxl)
library(ggrepel)

#### Inserting functions ####

qcut2 <- function(x, n) {
  findInterval(x, quantile(x, seq(0, 1, length = n + 1), na.rm=T), all.inside = T)
}


#### Loading DATASETS ####
issp_2007 <- read_dta("Datasets/ZA4850_v2-0-0.dta")
issp_2007$V5 <- as.numeric(issp_2007$V5)

#### Checking and Preparing Sleep Variables in the DATASETS ####
table(issp_2007$V67, exclude = NULL)
table(issp_2007$V68, exclude = NULL)

issp_2007$sleeping_time_hour_alternative <- substr(as.POSIXct(sprintf("%04.0f", issp_2007$V68), format='%H%M'), 12, 16)
issp_2007$getting_up_time_hour_alternative <- substr(as.POSIXct(sprintf("%04.0f", issp_2007$V67), format='%H%M'), 12, 16)

table(issp_2007$sleeping_time_hour_alternative, exclude = NULL)
table(issp_2007$getting_up_time_hour_alternative, exclude = NULL)

issp_2007$Sleeping_Hour <- format(as.POSIXct(issp_2007$sleeping_time_hour_alternative, format="%H:%M"),"%H")
issp_2007$Sleeping_Minute <- format(as.POSIXct(issp_2007$sleeping_time_hour_alternative, format="%H:%M"),"%M")

issp_2007$Waking_Hour <- format(as.POSIXct(issp_2007$getting_up_time_hour_alternative, format="%H:%M"),"%H")
issp_2007$Waking_Minute <- format(as.POSIXct(issp_2007$getting_up_time_hour_alternative, format="%H:%M"),"%M")

issp_2007$sleeping_time_hour_alternative <- as.numeric(issp_2007$Sleeping_Hour) + as.numeric(issp_2007$Sleeping_Minute)/60
issp_2007$getting_up_time_hour_alternative <- as.numeric(issp_2007$Waking_Hour) + as.numeric(issp_2007$Waking_Minute)/60

#### Loading Netherlands Dataset ####
issp_2007_netherlands <- read_dta("Datasets/issp_2007nl_def.dta")

issp_2007_netherlands$nl_hinc <- ifelse(issp_2007_netherlands$nl_hinc==99999, NA, issp_2007_netherlands$nl_hinc)
issp_2007_netherlands$NL_INC <- issp_2007_netherlands$nl_hinc
table(issp_2007_netherlands$NL_INC, exclude = NULL)
str(issp_2007_netherlands$NL_INC)
issp_2007_netherlands$NL_INC <- as.numeric(issp_2007_netherlands$NL_INC)
table(qcut2(issp_2007_netherlands$NL_INC, 5), exclude = NULL)
issp_2007_netherlands$NL_INC_rank <- qcut2(issp_2007_netherlands$NL_INC, 5)
table(issp_2007_netherlands$NL_INC_rank, exclude = NULL)

issp_2007_netherlands$V58 <- ifelse(issp_2007_netherlands$V58==-1 | issp_2007_netherlands$V58==8, NA, issp_2007_netherlands$V58)
table(issp_2007_netherlands$V58, exclude = NULL)

issp_2007_netherlands$V66 <- ifelse(issp_2007_netherlands$V66==-1, NA, issp_2007_netherlands$V66)
table(issp_2007_netherlands$V66, exclude = NULL)

issp_2007_netherlands$V67 <- ifelse(issp_2007_netherlands$V67==-1, NA, issp_2007_netherlands$V67)
table(issp_2007_netherlands$V67, exclude = NULL)

issp_2007_netherlands$V68 <- ifelse(issp_2007_netherlands$V68==-1, NA, issp_2007_netherlands$V68)
table(issp_2007_netherlands$V68, exclude = NULL)

issp_2007_netherlands$degree <- ifelse(issp_2007_netherlands$degree==9, NA, issp_2007_netherlands$degree)
table(issp_2007_netherlands$degree, exclude = NULL)

issp_2007_netherlands$nl_hinc <- ifelse(issp_2007_netherlands$nl_hinc==99999, NA, issp_2007_netherlands$nl_hinc)
table(issp_2007_netherlands$nl_hinc, exclude = NULL)

issp_2007_netherlands$PARTY_LR <- ifelse(issp_2007_netherlands$PARTY_LR==7 | issp_2007_netherlands$PARTY_LR==9, NA, issp_2007_netherlands$PARTY_LR)
table(issp_2007_netherlands$PARTY_LR, exclude = NULL)

issp_2007_netherlands$attend <- ifelse(issp_2007_netherlands$attend==98 | issp_2007_netherlands$attend==99, NA, issp_2007_netherlands$attend)
table(issp_2007_netherlands$attend, exclude = NULL)

table(issp_2007_netherlands$relig, exclude = NULL)
issp_2007_netherlands$relig <- ifelse(issp_2007_netherlands$relig==0, NA, issp_2007_netherlands$relig)

table(issp_2007_netherlands$religgrp, exclude = NULL)
issp_2007_netherlands$religgrp <- ifelse(issp_2007_netherlands$religgrp==99, NA, issp_2007_netherlands$religgrp)

#### Checking and Preparing Sleep Variables in the Netherlands Dataset ####

issp_2007_netherlands$V67 <- as.numeric(issp_2007_netherlands$V67)
table(issp_2007_netherlands$V67, exclude = NULL)
issp_2007_netherlands$V67 <- issp_2007_netherlands$V67*100

issp_2007_netherlands$V68 <- as.numeric(issp_2007_netherlands$V68)
table(issp_2007_netherlands$V68, exclude = NULL)
issp_2007_netherlands$V68 <- issp_2007_netherlands$V68*100

issp_2007_netherlands$sleeping_time_hour_alternative <- substr(as.POSIXct(sprintf("%04.0f", issp_2007_netherlands$V68), format='%H%M'), 12, 16)
issp_2007_netherlands$getting_up_time_hour_alternative <- substr(as.POSIXct(sprintf("%04.0f", issp_2007_netherlands$V67), format='%H%M'), 12, 16)

table(issp_2007_netherlands$sleeping_time_hour_alternative, exclude = NULL)
table(issp_2007_netherlands$getting_up_time_hour_alternative, exclude = NULL)

issp_2007_netherlands$Sleeping_Hour <- format(as.POSIXct(issp_2007_netherlands$sleeping_time_hour_alternative, format="%H:%M"),"%H")
issp_2007_netherlands$Sleeping_Minute <- format(as.POSIXct(issp_2007_netherlands$sleeping_time_hour_alternative, format="%H:%M"),"%M")

issp_2007_netherlands$Waking_Hour <- format(as.POSIXct(issp_2007_netherlands$getting_up_time_hour_alternative, format="%H:%M"),"%H")
issp_2007_netherlands$Waking_Minute <- format(as.POSIXct(issp_2007_netherlands$getting_up_time_hour_alternative, format="%H:%M"),"%M")

issp_2007_netherlands$sleeping_time_hour_alternative <- as.numeric(issp_2007_netherlands$Sleeping_Hour) + as.numeric(issp_2007_netherlands$Sleeping_Minute)/60
issp_2007_netherlands$getting_up_time_hour_alternative <- as.numeric(issp_2007_netherlands$Waking_Hour) + as.numeric(issp_2007_netherlands$Waking_Minute)/60

#### Merging ISSP 2007 and the Netherlands Datasets ####
issp_2007 <- dplyr::bind_rows(issp_2007, issp_2007_netherlands)

table(issp_2007$religgrp, exclude = NULL)

issp_2007$religious_denomination <- NA
issp_2007$religious_denomination[issp_2007$religgrp==1] <- "No_religion"
issp_2007$religious_denomination[issp_2007$religgrp==2] <- "Roman_Catholic"
issp_2007$religious_denomination[issp_2007$religgrp==3] <- "Protestant"
issp_2007$religious_denomination[issp_2007$religgrp==4] <- "Christian_Orthodox"
issp_2007$religious_denomination[issp_2007$religgrp==5] <- "Jewish"
issp_2007$religious_denomination[issp_2007$religgrp==6] <- "Islam"
issp_2007$religious_denomination[issp_2007$religgrp==7] <- "Buddhism"
issp_2007$religious_denomination[issp_2007$religgrp==8] <- "Hinduism"
issp_2007$religious_denomination[issp_2007$religgrp==9] <- "Other_Christian_Religions"
issp_2007$religious_denomination[issp_2007$religgrp==10] <- "Other_Eastern_Religions"
issp_2007$religious_denomination[issp_2007$religgrp==11] <- "Other_Religions"
table(issp_2007$religious_denomination, exclude = NULL)

issp_2007 <- fastDummies::dummy_cols(issp_2007, select_columns = 'religious_denomination')

table(issp_2007$religious_denomination_Buddhism, exclude = NULL)
issp_2007$religious_denomination_Buddhism_avg <- mean(issp_2007$religious_denomination_Buddhism, na.rm = T)
issp_2007$religious_denomination_Buddhism_beta <- (issp_2007$religious_denomination_Buddhism-issp_2007$religious_denomination_Buddhism_avg)
table(issp_2007$religious_denomination_Buddhism_beta, exclude = NULL)

table(issp_2007$religious_denomination_Christian_Orthodox, exclude = NULL)
issp_2007$religious_denomination_Christian_Orthodox_avg <- mean(issp_2007$religious_denomination_Christian_Orthodox, na.rm = T)
issp_2007$religious_denomination_Christian_Orthodox_beta <- (issp_2007$religious_denomination_Christian_Orthodox-issp_2007$religious_denomination_Christian_Orthodox_avg)
table(issp_2007$religious_denomination_Christian_Orthodox_beta, exclude = NULL)

table(issp_2007$religious_denomination_Hinduism, exclude = NULL)
issp_2007$religious_denomination_Hinduism_avg <- mean(issp_2007$religious_denomination_Hinduism, na.rm = T)
issp_2007$religious_denomination_Hinduism_beta <- (issp_2007$religious_denomination_Hinduism-issp_2007$religious_denomination_Hinduism_avg)
table(issp_2007$religious_denomination_Hinduism_beta, exclude = NULL)

table(issp_2007$religious_denomination_Islam, exclude = NULL)
issp_2007$religious_denomination_Islam_avg <- mean(issp_2007$religious_denomination_Islam, na.rm = T)
issp_2007$religious_denomination_Islam_beta <- (issp_2007$religious_denomination_Islam-issp_2007$religious_denomination_Islam_avg)
table(issp_2007$religious_denomination_Islam_beta, exclude = NULL)

table(issp_2007$religious_denomination_Jewish, exclude = NULL)
issp_2007$religious_denomination_Jewish_avg <- mean(issp_2007$religious_denomination_Jewish, na.rm = T)
issp_2007$religious_denomination_Jewish_beta <- (issp_2007$religious_denomination_Jewish-issp_2007$religious_denomination_Jewish_avg)
table(issp_2007$religious_denomination_Jewish_beta, exclude = NULL)

table(issp_2007$religious_denomination_No_religion, exclude = NULL)
issp_2007$religious_denomination_No_religion_avg <- mean(issp_2007$religious_denomination_No_religion, na.rm = T)
issp_2007$religious_denomination_No_religion_beta <- (issp_2007$religious_denomination_No_religion-issp_2007$religious_denomination_No_religion_avg)
table(issp_2007$religious_denomination_No_religion_beta, exclude = NULL)

table(issp_2007$religious_denomination_Other_Christian_Religions, exclude = NULL)
issp_2007$religious_denomination_Other_Christian_Religions_avg <- mean(issp_2007$religious_denomination_Other_Christian_Religions, na.rm = T)
issp_2007$religious_denomination_Other_Christian_Religions_beta <- (issp_2007$religious_denomination_Other_Christian_Religions-issp_2007$religious_denomination_Other_Christian_Religions_avg)
table(issp_2007$religious_denomination_Other_Christian_Religions_beta, exclude = NULL)

table(issp_2007$religious_denomination_Other_Eastern_Religions, exclude = NULL)
issp_2007$religious_denomination_Other_Eastern_Religions_avg <- mean(issp_2007$religious_denomination_Other_Eastern_Religions, na.rm = T)
issp_2007$religious_denomination_Other_Eastern_Religions_beta <- (issp_2007$religious_denomination_Other_Eastern_Religions-issp_2007$religious_denomination_Other_Eastern_Religions_avg)
table(issp_2007$religious_denomination_Other_Eastern_Religions_beta, exclude = NULL)

table(issp_2007$religious_denomination_Other_Religions, exclude = NULL)
issp_2007$religious_denomination_Other_Religions_avg <- mean(issp_2007$religious_denomination_Other_Religions, na.rm = T)
issp_2007$religious_denomination_Other_Religions_beta <- (issp_2007$religious_denomination_Other_Religions-issp_2007$religious_denomination_Other_Religions_avg)
table(issp_2007$religious_denomination_Other_Religions_beta, exclude = NULL)

table(issp_2007$religious_denomination_Protestant, exclude = NULL)
issp_2007$religious_denomination_Protestant_avg <- mean(issp_2007$religious_denomination_Protestant, na.rm = T)
issp_2007$religious_denomination_Protestant_beta <- (issp_2007$religious_denomination_Protestant-issp_2007$religious_denomination_Protestant_avg)
table(issp_2007$religious_denomination_Protestant_beta, exclude = NULL)

table(issp_2007$religious_denomination_Roman_Catholic, exclude = NULL)
issp_2007$religious_denomination_Roman_Catholic_avg <- mean(issp_2007$religious_denomination_Roman_Catholic, na.rm = T)
issp_2007$religious_denomination_Roman_Catholic_beta <- (issp_2007$religious_denomination_Roman_Catholic-issp_2007$religious_denomination_Roman_Catholic_avg)
table(issp_2007$religious_denomination_Roman_Catholic_beta, exclude = NULL)

#### Dropping Observations with Missing Sleep Variables ####
issp_2007 <- subset(issp_2007, subset = (!is.na(issp_2007$V67) | !is.na(issp_2007$V68)))

#### Calculating Duration of Sleep ####
issp_2007$sleep_duration_alternative <-  (issp_2007$getting_up_time_hour_alternative-issp_2007$sleeping_time_hour_alternative) 
table(issp_2007$sleep_duration_alternative, exclude = NULL)
issp_2007$sleep_duration_alternative <-  ifelse(issp_2007$getting_up_time_hour_alternative>=issp_2007$sleeping_time_hour_alternative, issp_2007$sleep_duration_alternative, (24 - issp_2007$sleeping_time_hour_alternative) + issp_2007$getting_up_time_hour_alternative) 
table(issp_2007$sleep_duration_alternative, exclude = NULL)

#### Calculating Chronotype ####
issp_2007$mid_point_alternative <- issp_2007$getting_up_time_hour_alternative - issp_2007$sleep_duration_alternative/2
issp_2007$mid_point_alternative <- ifelse(issp_2007$mid_point_alternative<0, issp_2007$mid_point_alternative+24, issp_2007$mid_point_alternative)
table(issp_2007$mid_point_alternative)

#### Exclusion Criteria based on Duration of Sleep and Chronotype ####
#issp_2007 <- subset(issp_2007, subset = (issp_2007$sleeping_time_hour_alternative<=6 | issp_2007$sleeping_time_hour_alternative>=19) & (issp_2007$getting_up_time_hour_alternative<=15 & issp_2007$getting_up_time_hour_alternative>=4) & (issp_2007$sleep_duration_alternative>=2 & sleep_duration_alternative<=18))

issp_2007$EXCLUDED <- NA
issp_2007$EXCLUDED <- ifelse((issp_2007$sleeping_time_hour_alternative<=6 | issp_2007$sleeping_time_hour_alternative>=19) & (issp_2007$getting_up_time_hour_alternative<=15 & issp_2007$getting_up_time_hour_alternative>=4) & (issp_2007$sleep_duration_alternative>=2 & issp_2007$sleep_duration_alternative<=18), "Included", "Excluded")
table(issp_2007$EXCLUDED, exclude = NULL)

#### Reverse Coding Chronotype - Bigger Values for Morningness and Smaller Values for Eveningness ####
table(issp_2007$mid_point_alternative, exclude = NULL)
issp_2007$ORIGINAL_MID_POINT <- issp_2007$mid_point_alternative
issp_2007$midcenter <- issp_2007$mid_point_alternative-16
table(issp_2007$midcenter, exclude = NULL)
issp_2007$cond3 <- issp_2007$midcenter + 24
issp_2007$midcenter <- ifelse(issp_2007$midcenter<0, issp_2007$cond3, issp_2007$midcenter)
table(issp_2007$midcenter, exclude = NULL)
issp_2007$mid_point_alternative <- issp_2007$midcenter
table(issp_2007$mid_point_alternative, exclude = NULL)
issp_2007$chronotype_interval <-  issp_2007$mid_point_alternative*-1
issp_2007$sleep_duration <-  issp_2007$sleep_duration_alternative
table(issp_2007$chronotype_interval, exclude = NULL)
table(issp_2007$sleep_duration, exclude = NULL)

#### Checking and Prepating Covariates ####
table(issp_2007$degree, exclude = NULL)
issp_2007$degree <- ifelse(issp_2007$degree>5, NA, issp_2007$degree)

table(issp_2007$sex, exclude = NULL) # 1 - male, 2 - female
issp_2007$sex <- issp_2007$sex
issp_2007$sex[issp_2007$sex==1] <- 1
issp_2007$sex[issp_2007$sex==2] <- 0
table(issp_2007$sex, exclude = NULL)


table(issp_2007$V66, exclude = NULL)
issp_2007$day_off <- NA
issp_2007$day_off[issp_2007$V66==2] <- 1 ## Weekend
issp_2007$day_off[issp_2007$V66==1] <- 0 ## Weekday
table(issp_2007$day_off, exclude = NULL)


table(issp_2007$V58, exclude = NULL)
issp_2007$political_interest <- issp_2007$V58
table(issp_2007$political_interest, exclude = NULL)
issp_2007$political_interest <- ifelse(issp_2007$political_interest>4, NA, issp_2007$political_interest)
issp_2007$political_interest <- 5 - issp_2007$political_interest

table(issp_2007$urbrural, exclude = NULL)
issp_2007$urban_rural_area <- NA
issp_2007$urban_rural_area[issp_2007$urbrural==5 | issp_2007$urbrural==4] <- 1 # to match with Greece data
issp_2007$urban_rural_area[issp_2007$urbrural==3] <- 2
issp_2007$urban_rural_area[issp_2007$urbrural==2] <- 3
issp_2007$urban_rural_area[issp_2007$urbrural==1] <- 4
table(issp_2007$urban_rural_area, exclude = NULL)

table(issp_2007$attend, exclude = NULL)

issp_2007$religiosity <- NA
issp_2007$religiosity[issp_2007$attend==8] <- 1
issp_2007$religiosity[issp_2007$attend==7 | issp_2007$attend==6] <- 2
issp_2007$religiosity[issp_2007$attend==5] <- 3
issp_2007$religiosity[issp_2007$attend==4] <- 4
issp_2007$religiosity[issp_2007$attend==3] <- 5
issp_2007$religiosity[issp_2007$attend==1 | issp_2007$attend==2] <- 6
table(issp_2007$religiosity, exclude = NULL)

issp_2007$age <- ifelse(issp_2007$age<18, NA, issp_2007$age)

issp_2007$PARTY_LR <- ifelse(issp_2007$PARTY_LR>5, NA, issp_2007$PARTY_LR)
table(issp_2007$PARTY_LR, exclude = NULL)

table(issp_2007$VOTE_LE, exclude = NULL)
issp_2007$turnout <- NA
issp_2007$turnout[issp_2007$VOTE_LE==2] <- 0
issp_2007$turnout[issp_2007$VOTE_LE==1] <- 1
table(issp_2007$turnout, exclude = NULL)

#### Loading Greece Dataset ####

greece_data <- read_dta("Datasets/Greece_2020.dta")

greece_data$V67 <-  greece_data$hourtime1 + greece_data$minutetime1/100
table(greece_data$V67, exclude = NULL)
greece_data$V67 <- greece_data$V67*100

greece_data$V68 <-  greece_data$hourtime2_corrected + greece_data$minutetime2/100
table(greece_data$V68, exclude = NULL)
greece_data$V68 <- greece_data$V68*100

greece_data$sleeping_time_hour_alternative <- substr(as.POSIXct(sprintf("%04.0f", greece_data$V68), format='%H%M'), 12, 16)
greece_data$getting_up_time_hour_alternative <- substr(as.POSIXct(sprintf("%04.0f", greece_data$V67), format='%H%M'), 12, 16)

table(greece_data$sleeping_time_hour_alternative, exclude = NULL)
table(greece_data$getting_up_time_hour_alternative, exclude = NULL)

greece_data$Sleeping_Hour <- format(as.POSIXct(greece_data$sleeping_time_hour_alternative, format="%H:%M"),"%H")
greece_data$Sleeping_Minute <- format(as.POSIXct(greece_data$sleeping_time_hour_alternative, format="%H:%M"),"%M")

greece_data$Waking_Hour <- format(as.POSIXct(greece_data$getting_up_time_hour_alternative, format="%H:%M"),"%H")
greece_data$Waking_Minute <- format(as.POSIXct(greece_data$getting_up_time_hour_alternative, format="%H:%M"),"%M")

greece_data$sleeping_time_hour_alternative <- as.numeric(greece_data$Sleeping_Hour) + as.numeric(greece_data$Sleeping_Minute)/60
greece_data$getting_up_time_hour_alternative <- as.numeric(greece_data$Waking_Hour) + as.numeric(greece_data$Waking_Minute)/60

greece_data$sleep_duration_alternative <-  (greece_data$getting_up_time_hour_alternative-greece_data$sleeping_time_hour_alternative) 
table(greece_data$sleep_duration_alternative, exclude = NULL)
greece_data$sleep_duration_alternative <-  ifelse(greece_data$getting_up_time_hour_alternative>=greece_data$sleeping_time_hour_alternative, greece_data$sleep_duration_alternative, (24 - greece_data$sleeping_time_hour_alternative) + greece_data$getting_up_time_hour_alternative) 
table(greece_data$sleep_duration_alternative, exclude = NULL)

greece_data$mid_point_alternative <- greece_data$getting_up_time_hour_alternative - greece_data$sleep_duration_alternative/2
greece_data$mid_point_alternative <- ifelse(greece_data$mid_point_alternative<0, greece_data$mid_point_alternative+24, greece_data$mid_point_alternative)
table(greece_data$mid_point_alternative)

#greece_data <- subset(greece_data, subset = (greece_data$sleeping_time_hour_alternative<=6 | greece_data$sleeping_time_hour_alternative>=19) & (greece_data$getting_up_time_hour_alternative<=15 & greece_data$getting_up_time_hour_alternative>=4) & (greece_data$sleep_duration_alternative>=2 & sleep_duration_alternative<=18))
greece_data$EXCLUDED <- NA
greece_data$EXCLUDED <- ifelse((greece_data$sleeping_time_hour_alternative<=6 | greece_data$sleeping_time_hour_alternative>=19) & (greece_data$getting_up_time_hour_alternative<=15 & greece_data$getting_up_time_hour_alternative>=4) & (greece_data$sleep_duration_alternative>=2 & greece_data$sleep_duration_alternative<=18), "Included", "Excluded")
table(greece_data$EXCLUDED, exclude = NULL)

table(greece_data$mid_point_alternative, exclude = NULL)
greece_data$ORIGINAL_MID_POINT <- greece_data$mid_point_alternative
greece_data$midcenter <- greece_data$mid_point_alternative-16
table(greece_data$midcenter, exclude = NULL)
greece_data$cond3 <- greece_data$midcenter + 24
greece_data$midcenter <- ifelse(greece_data$midcenter<0, greece_data$cond3, greece_data$midcenter)
table(greece_data$midcenter, exclude = NULL)
greece_data$mid_point_alternative <- greece_data$midcenter
table(greece_data$mid_point_alternative, exclude = NULL)
greece_data$chronotype_interval <-  greece_data$mid_point_alternative*-1
greece_data$sleep_duration <-  greece_data$sleep_duration_alternative
table(greece_data$chronotype_interval, exclude = NULL)
table(greece_data$sleep_duration, exclude = NULL)

greece_data$degree <- NA
greece_data$degree[greece_data$D03f==1] <- 0
greece_data$degree[greece_data$D03f==2 | greece_data$D03f==3] <- 1
greece_data$degree[greece_data$D03f==4 | greece_data$D03f==5] <- 2
greece_data$degree[greece_data$D03f==6] <- 3
greece_data$degree[greece_data$D03f==7] <- 4
greece_data$degree[greece_data$D03f==8 | greece_data$D03f==9 | greece_data$D03f==10] <- 5
table(greece_data$degree, exclude = NULL)

greece_data$sex <- NA
greece_data$sex[greece_data$D02f==1] <- 1
greece_data$sex[greece_data$D02f==2] <- 0
table(greece_data$sex, exclude = NULL)

greece_data$interview_time <- paste(greece_data$qyear, greece_data$qmonth, greece_data$qday, sep = "-")
table(greece_data$interview_time, exclude = NULL)

greece_data$day_off <- NA
greece_data$day_off[greece_data$interview_time=="2019-12-14" |
                      greece_data$interview_time=="2019-12-15" |
                      greece_data$interview_time== "2019-12-22" |
                      greece_data$interview_time=="2019-12-28" |
                      greece_data$interview_time=="2019-12-29" |
                      greece_data$interview_time== "2020-1-4" |
                      greece_data$interview_time=="2020-1-5" |
                      greece_data$interview_time== "2020-2-1" |
                      greece_data$interview_time== "2020-2-2" |
                      greece_data$interview_time=="2020-2-8" |
                      greece_data$interview_time=="2020-2-9" |
                      greece_data$interview_time=="2020-2-15" |
                      greece_data$interview_time== "2020-2-16"] <- 1
greece_data$day_off <- ifelse(is.na(greece_data$day_off), 0, greece_data$day_off)
table(greece_data$day_off, exclude = NULL)

table(greece_data$Q01f, exclude = NULL)
greece_data$political_interest <- greece_data$Q01f
table(greece_data$political_interest, exclude = NULL)
greece_data$political_interest <- ifelse(greece_data$political_interest>4, NA, greece_data$political_interest)
greece_data$political_interest <- 5 - greece_data$political_interest

greece_data$urban_rural_area <- greece_data$D19f
greece_data$urban_rural_area <- ifelse(greece_data$urban_rural_area==5, NA, greece_data$urban_rural_area)
table(greece_data$urban_rural_area, exclude = NULL)

table(greece_data$D11f, exclude = NULL)
greece_data$religiosity <- greece_data$D11f
greece_data$religiosity <- ifelse(greece_data$religiosity==7, NA, greece_data$religiosity)
table(greece_data$religiosity, exclude = NULL)

table(greece_data$Q18f, exclude = NULL)
greece_data$ideology <- greece_data$Q18f
greece_data$ideology <- ifelse(greece_data$ideology==12, NA, greece_data$ideology)
table(greece_data$ideology, exclude = NULL)
greece_data$PARTY_LR <- (5-1)/(11-1)*(greece_data$ideology-11) + 5
table(greece_data$PARTY_LR, exclude = NULL)
greece_data$PARTY_LR <- round(greece_data$PARTY_LR, 0)
greece_data$PARTY_LR[greece_data$ideology==5] <- 2
greece_data$PARTY_LR[greece_data$ideology==7] <- 4
table(greece_data$PARTY_LR, greece_data$ideology, exclude = NULL)
table(greece_data$PARTY_LR, exclude = NULL)

table(greece_data$age, exclude = NULL)
greece_data$age <- ifelse(greece_data$age<18, NA, greece_data$age)
greece_data <- subset(greece_data, subset = !is.na(age))

table(greece_data$Q12LHaf, exclude = NULL)
greece_data$turnout <- NA
greece_data$turnout[greece_data$Q12LHaf==2] <- 0
greece_data$turnout[greece_data$Q12LHaf==1] <- 1
table(greece_data$turnout, exclude = NULL)

table(greece_data$D09f, exclude = NULL)
str(greece_data$D09f)
greece_data$D09f <- as.numeric(greece_data$D09f)
greece_data$GR_INC_rank <- greece_data$D09f
greece_data$GR_INC_rank <- ifelse(greece_data$GR_INC_rank==6, NA, greece_data$GR_INC_rank)
table(greece_data$GR_INC_rank, exclude = NULL)

greece_data$V5 <- 2020

issp_2007 <- dplyr::bind_rows(issp_2007, greece_data)

#### Income variable for all ####

table(issp_2007$FI_INC, exclude = NULL)
str(issp_2007$FI_INC)
issp_2007$FI_INC <- as.numeric(issp_2007$FI_INC)
table(qcut2(issp_2007$FI_INC, 5), exclude = NULL)
issp_2007$FI_INC_rank <- qcut2(issp_2007$FI_INC, 5)
table(issp_2007$FI_INC_rank, exclude = NULL)

table(issp_2007$IE_INC, exclude = NULL)
str(issp_2007$IE_INC)
issp_2007$IE_INC <- as.numeric(issp_2007$IE_INC)
table(qcut2(issp_2007$IE_INC, 5), exclude = NULL)
issp_2007$IE_INC_rank <- qcut2(issp_2007$IE_INC, 5)
table(issp_2007$IE_INC_rank, exclude = NULL)

table(issp_2007$KR_INC, exclude = NULL)
str(issp_2007$KR_INC)
issp_2007$KR_INC <- as.numeric(issp_2007$KR_INC)
table(qcut2(issp_2007$KR_INC, 5), exclude = NULL)
issp_2007$KR_INC_rank <- qcut2(issp_2007$KR_INC, 5)
table(issp_2007$KR_INC_rank, exclude = NULL)

table(issp_2007$MX_INC, exclude = NULL)
str(issp_2007$MX_INC)
issp_2007$MX_INC <- as.numeric(issp_2007$MX_INC)
table(qcut2(issp_2007$MX_INC, 5), exclude = NULL)
issp_2007$MX_INC_rank <- qcut2(issp_2007$MX_INC, 5)
table(issp_2007$MX_INC_rank, exclude = NULL)

table(issp_2007$PH_INC, exclude = NULL)
str(issp_2007$PH_INC)
issp_2007$PH_INC <- as.numeric(issp_2007$PH_INC)
table(qcut2(issp_2007$PH_INC, 5), exclude = NULL)
issp_2007$PH_INC_rank <- qcut2(issp_2007$PH_INC, 5)
table(issp_2007$PH_INC_rank, exclude = NULL)

table(issp_2007$RU_INC, exclude = NULL)
str(issp_2007$RU_INC)
issp_2007$RU_INC <- as.numeric(issp_2007$RU_INC)
table(qcut2(issp_2007$RU_INC, 5), exclude = NULL)
issp_2007$RU_INC_rank <- qcut2(issp_2007$RU_INC, 5)
table(issp_2007$RU_INC_rank, exclude = NULL)

table(issp_2007$NZ_INC, exclude = NULL)
str(issp_2007$NZ_INC)
issp_2007$NZ_INC <- as.numeric(issp_2007$NZ_INC)
table(qcut2(issp_2007$NZ_INC, 5), exclude = NULL)
issp_2007$NZ_INC_rank <- qcut2(issp_2007$NZ_INC, 5)
table(issp_2007$NZ_INC_rank, exclude = NULL)

issp_2007$income <- rowSums(issp_2007[,c("FI_INC_rank", 
                                         "IE_INC_rank",
                                         "KR_INC_rank",
                                         "MX_INC_rank",
                                         "NL_INC_rank",
                                         "NZ_INC_rank",
                                         "PH_INC_rank",
                                         "RU_INC_rank",
                                         "GR_INC_rank")], na.rm = T)
table(issp_2007$income, exclude = NULL)
table(issp_2007$income, issp_2007$FI_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$IE_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$KR_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$MX_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$NL_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$NZ_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$PH_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$RU_INC_rank, exclude = NULL)
table(issp_2007$income, issp_2007$GR_INC_rank, exclude = NULL)
issp_2007$income <- ifelse(issp_2007$income==0, NA, issp_2007$income)
table(issp_2007$income, exclude = NULL)

table(issp_2007$V5, exclude = NULL)
issp_2007$country <- NA
issp_2007$country[issp_2007$V5==246] <- 'Finland'
issp_2007$country[issp_2007$V5==2020] <- 'Greece'
issp_2007$country[issp_2007$V5==372] <- 'Ireland'
issp_2007$country[issp_2007$V5==410] <- 'South Korea'
issp_2007$country[issp_2007$V5==484] <- 'Mexico'
issp_2007$country[issp_2007$V5==528] <- 'the Netherlands'
issp_2007$country[issp_2007$V5==554] <- 'New Zealand'
issp_2007$country[issp_2007$V5==608] <- 'the Philippines'
issp_2007$country[issp_2007$V5==643] <- 'Russia'
table(issp_2007$country, exclude = NULL)
issp_2007$country <- factor(issp_2007$country, 
                            levels = c('Finland','Greece', 'Ireland', 'Mexico', 'the Netherlands', 'New Zealand', 'the Philippines', 'Russia', 'South Korea'))
table(issp_2007$country, exclude = NULL)


issp_2007 <- subset(issp_2007, subset = country == 'Finland' | 
                      country == 'Greece' |
                      country == 'Ireland' | 
                      country == 'South Korea' | 
                      country == 'Mexico' | 
                      country == 'the Netherlands' | 
                      country == 'New Zealand' | 
                      country == 'the Philippines' | country == 'Russia')

issp_2007$local_election <- NA
issp_2007$local_election[issp_2007$country=="the Netherlands" | issp_2007$country=="South Korea"] <- 1
issp_2007$local_election <- ifelse(is.na(issp_2007$local_election), 0, issp_2007$local_election)
table(issp_2007$local_election, exclude = NULL)

issp_2007$polyarchy <- NA
issp_2007$polyarchy[issp_2007$country=="Finland"] <- .887
issp_2007$polyarchy[issp_2007$country=="Greece"] <- .862
issp_2007$polyarchy[issp_2007$country=="Ireland"] <- .888
issp_2007$polyarchy[issp_2007$country=="Mexico"] <- .672
issp_2007$polyarchy[issp_2007$country=="the Netherlands"] <- .864
issp_2007$polyarchy[issp_2007$country=="New Zealand"] <- .887
issp_2007$polyarchy[issp_2007$country=="the Philippines"] <- .498
issp_2007$polyarchy[issp_2007$country=="Russia"] <- .324
issp_2007$polyarchy[issp_2007$country=="South Korea"] <- .85
table(issp_2007$polyarchy, exclude = NULL)

#### Getting Country-Region Level Variables ####
table(issp_2007$IE_REG, exclude = NULL)
table(issp_2007$FI_REG, exclude = NULL)
table(issp_2007$KR_REG, exclude = NULL)
table(issp_2007$MX_REG, exclude = NULL)
table(issp_2007$NL_REG, exclude = NULL)
table(issp_2007$NZ_REG, exclude = NULL)
table(issp_2007$PH_REG, exclude = NULL)
table(issp_2007$RU_REG, exclude = NULL)

issp_2007$region_id <- NA
issp_2007$region_id <- rowSums( issp_2007[,c("IE_REG", "FI_REG", "KR_REG",
                                             "MX_REG", "NL_REG", "NZ_REG",
                                             "PH_REG", "RU_REG")], na.rm = T)

table(issp_2007$region_id, exclude = NULL)

issp_2007$region_id <- ifelse(issp_2007$region_id==0, NA, issp_2007$region_id)

issp_2007$region_id_merging <- paste(issp_2007$V5, issp_2007$region_id)
table(issp_2007$region_id_merging, exclude = NULL)  

geocodes_temperatures <- read_excel("Geographic Dataset/geocodes_temperatures.xlsx")
geocodes_temperatures$region_id_merging <- paste(geocodes_temperatures$V5, geocodes_temperatures$region_id)
table(geocodes_temperatures$region_id_merging, exclude = NULL) 

issp_2007 <- merge(issp_2007,geocodes_temperatures,by="region_id_merging", all = T)


#### Excluded and Included Sample Characteristics ####

#### Online Appendices, Fig. 1. Sleep patterns across the excluded and included samples.  ####

table(issp_2007$EXCLUDED, exclude = NULL)
issp_2007$EXCLUDED <- ifelse(is.na(issp_2007$EXCLUDED), "Missing Sleep-Wake Meaures", issp_2007$EXCLUDED)
issp_2007 <- subset(issp_2007, subset = EXCLUDED!="Missing Sleep-Wake Meaures")
table(issp_2007$EXCLUDED, exclude = NULL)

table(issp_2007$ORIGINAL_MID_POINT, exclude = NULL)
sort(unique(issp_2007$ORIGINAL_MID_POINT))
table(issp_2007$sleep_duration, exclude = NULL)

table(issp_2007$ORIGINAL_MID_POINT, exclude = NULL)
psych::describe(issp_2007$ORIGINAL_MID_POINT)

library(scales)

excluded <- subset(issp_2007, EXCLUDED =="Excluded", select = c("sleep_duration", "ORIGINAL_MID_POINT", "country", "getting_up_time_hour_alternative", "sleeping_time_hour_alternative", "sleep_duration_alternative", "turnout"))
included <- subset(issp_2007, EXCLUDED =="Included", select = c("sleep_duration", "ORIGINAL_MID_POINT", "country", "getting_up_time_hour_alternative", "sleeping_time_hour_alternative", "sleep_duration_alternative", "turnout"))

psych::describe(excluded$sleep_duration)

g.top.excluded <- ggplot(excluded, aes(x = sleep_duration, y = (..count..)/sum(..count..)*100)) +
  geom_histogram(color="black", alpha=0.6, position = 'identity', bins = 50) +
  scale_x_continuous(position = "top") +
  scale_y_reverse() +
  theme_classic() +
  theme(axis.text.x=element_text(colour="black", size = 10, angle = 90, vjust = 1, hjust=1), axis.text.y=element_text(colour="black", size = 10)) +
  theme(axis.title=element_text(size=10, colour="black", face="bold"), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
  labs(    x = "\nSleep Duration\n", 
           y = "\n Percentage (%) \n in the Whole Excluded Sample (N=622)\n",
           title = "") +
  theme(strip.text.x = element_text(size = 12, colour = "black")) +
  ggpubr::grids(linetype = "dashed") +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24), limits = c(0, 24),
                     labels = c("0 hr", "", "", "3 hrs", "", "", "6 hrs", "", "", "9 hrs", "", "", "12 hrs", "", "", "15 hrs", "", "", "18 hrs", "", "", "21 hrs", "", "", "24 hrs")) +
  scale_y_continuous(breaks = c(0, 2, 4, 6), limits = c(0, 6)) + facet_wrap(.~country) +
  theme(strip.text.x = element_text(size = 10, margin = margin (1.5, 1.5, 1.5, 1.5)))

g.top.excluded

psych::describe(excluded$ORIGINAL_MID_POINT)

g.bottom.excluded <- ggplot(excluded, aes(x = ORIGINAL_MID_POINT, y = (..count..)/sum(..count..)*100)) +
  geom_histogram(color="black", alpha=0.6, position = 'identity', bins = 50) +
  scale_x_continuous() +
  theme_classic() +
  theme(axis.text.x=element_text(colour="black", size = 10, angle = 90, vjust = 1, hjust=1), axis.text.y=element_text(colour="black", size = 10)) +
  theme(axis.title=element_text(size=10, colour="black", face="bold"), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
  labs(    x = "\nMid-Point Timing of Sleep\n", 
           y = "\n Percentage (%) \n in the Whole Excluded Sample (N=622)\n",
           title = "") +
  theme(strip.text.x = element_text(size = 12, colour = "black")) +
  ggpubr::grids(linetype = "dashed") +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 23.99), limits = c(0, 23.99),
                     labels = c("0:00", "", "", "3:00", "", "", "6:00", "", "", "9:00", "", "", "12:00", "", "", "15:00", "", "", "18:00", "", "", "21:00", "", "", "23:59")) +
  scale_y_continuous(breaks = c(0, 2, 4, 6), limits = c(0, 6)) + facet_wrap(.~country) +
  theme(strip.text.x = element_text(size = 10, margin = margin (1.5, 1.5, 1.5, 1.5)))

g.bottom.excluded

psych::describe(included$sleep_duration)

g.top.included <- ggplot(included, aes(x = sleep_duration, y = (..count..)/sum(..count..)*100)) +
  geom_histogram(color="black", alpha=0.6, position = 'identity', bins = 50) +
  theme_classic() +
  theme(axis.text.x=element_text(colour="black", size = 10, angle = 90, vjust = 1, hjust=1), axis.text.y=element_text(colour="black", size = 10)) +
  theme(axis.title=element_text(size=10, colour="black", face="bold"), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
  labs(    x = "\nSleep Duration\n", 
           y = "\n Percentage (%) \n in the Whole Included Sample (N=11,556)\n",
           title = "") +
  theme(strip.text.x = element_text(size = 12, colour = "black")) +
  ggpubr::grids(linetype = "dashed") +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24), limits = c(0, 24),
                     labels = c("0 hr", "", "", "3 hrs", "", "", "6 hrs", "", "", "9 hrs", "", "", "12 hrs", "", "", "15 hrs", "", "", "18 hrs", "", "", "21 hrs", "", "", "24 hrs")) +
  scale_y_continuous(breaks = c(0, 2, 4, 6), limits = c(0, 6)) + facet_wrap(.~country) +
  theme(strip.text.x = element_text(size = 10, margin = margin (1.5, 1.5, 1.5, 1.5)))

g.top.included

psych::describe(included$ORIGINAL_MID_POINT)

g.bottom.included <- ggplot(included, aes(x = ORIGINAL_MID_POINT, y = (..count..)/sum(..count..)*100)) +
  geom_histogram(color="black", alpha=0.6, position = 'identity', bins = 50) +
  scale_x_continuous() +
  theme_classic() +
  theme(axis.text.x=element_text(colour="black", size = 10, angle = 90, vjust = 1, hjust=1), axis.text.y=element_text(colour="black", size = 10)) +
  theme(axis.title=element_text(size=10, colour="black", face="bold"), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
  labs(    x = "\nMid-Point Timing of Sleep\n", 
           y = "\n Percentage (%) \n in the Whole Included Sample (N=11,556)\n",
           title = "") +
  theme(strip.text.x = element_text(size = 12, colour = "black")) +
  ggpubr::grids(linetype = "dashed") +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 23.99), limits = c(0, 23.99),
                     labels = c("0:00", "", "", "3:00", "", "", "6:00", "", "", "9:00", "", "", "12:00", "", "", "15:00", "", "", "18:00", "", "", "21:00", "", "", "23:59")) +
  scale_y_continuous(breaks = c(0, 2, 4, 6), limits = c(0, 6)) + facet_wrap(.~country) +
  theme(strip.text.x = element_text(size = 10, margin = margin (1.5, 1.5, 1.5, 1.5)))

g.bottom.included

library(cowplot)

excluded_included_sample_sleep <- plot_grid(g.top.excluded, g.top.included, g.bottom.excluded, g.bottom.included, 
                                            nrow = 2, align = 'hv', greedy = T)
ggsave(filename = "excluded_included_sample_sleep.jpg", plot = excluded_included_sample_sleep, width = 16, height = 8, dpi = 1000)

#### Online Appendices, Table 13 Three-way table summarizing the distribution of the exclusion criteria across the excluded sample’s sleep measures. ####

excluded$ST <- NA
excluded$ST[excluded$sleeping_time_hour_alternative>6 | excluded$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded$ST <- ifelse(is.na(excluded$ST), "Keep", excluded$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded$ST, exclude = NULL)


excluded$GT <- NA
excluded$GT[excluded$getting_up_time_hour_alternative>15 | excluded$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded$GT <- ifelse(is.na(excluded$GT), "Keep", excluded$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded$GT, exclude = NULL)


excluded$SD <- NA
excluded$SD[excluded$sleep_duration_alternative>18 | excluded$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded$SD <- ifelse(is.na(excluded$SD), "Keep", excluded$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded$SD, exclude = NULL)

excluded.table <- xtabs(~ ST + GT + SD, data=excluded)
excluded.table 
ftable(excluded.table)

excluded.finland <- subset(excluded, subset = country == "Finland")

excluded.finland$ST <- NA
excluded.finland$ST[excluded.finland$sleeping_time_hour_alternative>6 | excluded.finland$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.finland$ST <- ifelse(is.na(excluded.finland$ST), "Keep", excluded.finland$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.finland$ST, exclude = NULL)


excluded.finland$GT <- NA
excluded.finland$GT[excluded.finland$getting_up_time_hour_alternative>15 | excluded.finland$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.finland$GT <- ifelse(is.na(excluded.finland$GT), "Keep", excluded.finland$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.finland$GT, exclude = NULL)


excluded.finland$SD <- NA
excluded.finland$SD[excluded.finland$sleep_duration_alternative>18 | excluded.finland$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.finland$SD <- ifelse(is.na(excluded.finland$SD), "Keep", excluded.finland$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.finland$SD, exclude = NULL)

excluded.finland.table <- xtabs(~ ST + GT + SD, data=excluded.finland)
excluded.finland.table 
ftable(excluded.finland.table)

excluded.greece <- subset(excluded, subset = country == "Greece")

excluded.greece$ST <- NA
excluded.greece$ST[excluded.greece$sleeping_time_hour_alternative>6 | excluded.greece$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.greece$ST <- ifelse(is.na(excluded.greece$ST), "Keep", excluded.greece$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.greece$ST, exclude = NULL)


excluded.greece$GT <- NA
excluded.greece$GT[excluded.greece$getting_up_time_hour_alternative>15 | excluded.greece$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.greece$GT <- ifelse(is.na(excluded.greece$GT), "Keep", excluded.greece$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.greece$GT, exclude = NULL)


excluded.greece$SD <- NA
excluded.greece$SD[excluded.greece$sleep_duration_alternative>18 | excluded.greece$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.greece$SD <- ifelse(is.na(excluded.greece$SD), "Keep", excluded.greece$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.greece$SD, exclude = NULL)

excluded.greece.table <- xtabs(~ ST + GT + SD, data=excluded.greece)
excluded.greece.table 
ftable(excluded.greece.table)

excluded.ireland <- subset(excluded, subset = country == "Ireland")

excluded.ireland$ST <- NA
excluded.ireland$ST[excluded.ireland$sleeping_time_hour_alternative>6 | excluded.ireland$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.ireland$ST <- ifelse(is.na(excluded.ireland$ST), "Keep", excluded.ireland$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.ireland$ST, exclude = NULL)


excluded.ireland$GT <- NA
excluded.ireland$GT[excluded.ireland$getting_up_time_hour_alternative>15 | excluded.ireland$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.ireland$GT <- ifelse(is.na(excluded.ireland$GT), "Keep", excluded.ireland$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.ireland$GT, exclude = NULL)


excluded.ireland$SD <- NA
excluded.ireland$SD[excluded.ireland$sleep_duration_alternative>18 | excluded.ireland$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.ireland$SD <- ifelse(is.na(excluded.ireland$SD), "Keep", excluded.ireland$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.ireland$SD, exclude = NULL)

excluded.ireland.table <- xtabs(~ ST + GT + SD, data=excluded.ireland)
excluded.ireland.table 
ftable(excluded.ireland.table)

excluded.mexico <- subset(excluded, subset = country == "Mexico")

excluded.mexico$ST <- NA
excluded.mexico$ST[excluded.mexico$sleeping_time_hour_alternative>6 | excluded.mexico$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.mexico$ST <- ifelse(is.na(excluded.mexico$ST), "Keep", excluded.mexico$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.mexico$ST, exclude = NULL)


excluded.mexico$GT <- NA
excluded.mexico$GT[excluded.mexico$getting_up_time_hour_alternative>15 | excluded.mexico$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.mexico$GT <- ifelse(is.na(excluded.mexico$GT), "Keep", excluded.mexico$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.mexico$GT, exclude = NULL)


excluded.mexico$SD <- NA
excluded.mexico$SD[excluded.mexico$sleep_duration_alternative>18 | excluded.mexico$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.mexico$SD <- ifelse(is.na(excluded.mexico$SD), "Keep", excluded.mexico$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.mexico$SD, exclude = NULL)

excluded.mexico.table <- xtabs(~ ST + GT + SD, data=excluded.mexico)
excluded.mexico.table 
ftable(excluded.mexico.table)

excluded.netherlands <- subset(excluded, subset = country == "the Netherlands")

excluded.netherlands$ST <- NA
excluded.netherlands$ST[excluded.netherlands$sleeping_time_hour_alternative>6 | excluded.netherlands$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.netherlands$ST <- ifelse(is.na(excluded.netherlands$ST), "Keep", excluded.netherlands$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.netherlands$ST, exclude = NULL)


excluded.netherlands$GT <- NA
excluded.netherlands$GT[excluded.netherlands$getting_up_time_hour_alternative>15 | excluded.netherlands$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.netherlands$GT <- ifelse(is.na(excluded.netherlands$GT), "Keep", excluded.netherlands$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.netherlands$GT, exclude = NULL)


excluded.netherlands$SD <- NA
excluded.netherlands$SD[excluded.netherlands$sleep_duration_alternative>18 | excluded.netherlands$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.netherlands$SD <- ifelse(is.na(excluded.netherlands$SD), "Keep", excluded.netherlands$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.netherlands$SD, exclude = NULL)

excluded.netherlands.table <- xtabs(~ ST + GT + SD, data=excluded.netherlands)
excluded.netherlands.table 
ftable(excluded.netherlands.table)

excluded.new.zealand <- subset(excluded, subset = country == "New Zealand")

excluded.new.zealand$ST <- NA
excluded.new.zealand$ST[excluded.new.zealand$sleeping_time_hour_alternative>6 | excluded.new.zealand$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.new.zealand$ST <- ifelse(is.na(excluded.new.zealand$ST), "Keep", excluded.new.zealand$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.new.zealand$ST, exclude = NULL)


excluded.new.zealand$GT <- NA
excluded.new.zealand$GT[excluded.new.zealand$getting_up_time_hour_alternative>15 | excluded.new.zealand$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.new.zealand$GT <- ifelse(is.na(excluded.new.zealand$GT), "Keep", excluded.new.zealand$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.new.zealand$GT, exclude = NULL)


excluded.new.zealand$SD <- NA
excluded.new.zealand$SD[excluded.new.zealand$sleep_duration_alternative>18 | excluded.new.zealand$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.new.zealand$SD <- ifelse(is.na(excluded.new.zealand$SD), "Keep", excluded.new.zealand$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.new.zealand$SD, exclude = NULL)

excluded.new.zealand.table <- xtabs(~ ST + GT + SD, data=excluded.new.zealand)
excluded.new.zealand.table 
ftable(excluded.new.zealand.table)

excluded.philippines <- subset(excluded, subset = country == "the Philippines")

excluded.philippines$ST <- NA
excluded.philippines$ST[excluded.philippines$sleeping_time_hour_alternative>6 | excluded.philippines$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.philippines$ST <- ifelse(is.na(excluded.philippines$ST), "Keep", excluded.philippines$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.philippines$ST, exclude = NULL)


excluded.philippines$GT <- NA
excluded.philippines$GT[excluded.philippines$getting_up_time_hour_alternative>15 | excluded.philippines$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.philippines$GT <- ifelse(is.na(excluded.philippines$GT), "Keep", excluded.philippines$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.philippines$GT, exclude = NULL)


excluded.philippines$SD <- NA
excluded.philippines$SD[excluded.philippines$sleep_duration_alternative>18 | excluded.philippines$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.philippines$SD <- ifelse(is.na(excluded.philippines$SD), "Keep", excluded.philippines$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.philippines$SD, exclude = NULL)

excluded.philippines.table <- xtabs(~ ST + GT + SD, data=excluded.philippines)
excluded.philippines.table 
ftable(excluded.philippines.table)

excluded.russia <- subset(excluded, subset = country == "Russia")

excluded.russia$ST <- NA
excluded.russia$ST[excluded.russia$sleeping_time_hour_alternative>6 | excluded.russia$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.russia$ST <- ifelse(is.na(excluded.russia$ST), "Keep", excluded.russia$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.russia$ST, exclude = NULL)


excluded.russia$GT <- NA
excluded.russia$GT[excluded.russia$getting_up_time_hour_alternative>15 | excluded.russia$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.russia$GT <- ifelse(is.na(excluded.russia$GT), "Keep", excluded.russia$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.russia$GT, exclude = NULL)


excluded.russia$SD <- NA
excluded.russia$SD[excluded.russia$sleep_duration_alternative>18 | excluded.russia$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.russia$SD <- ifelse(is.na(excluded.russia$SD), "Keep", excluded.russia$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.russia$SD, exclude = NULL)

excluded.russia.table <- xtabs(~ ST + GT + SD, data=excluded.russia)
excluded.russia.table 
ftable(excluded.russia.table)

excluded.south.korea <- subset(excluded, subset = country == "South Korea")

excluded.south.korea$ST <- NA
excluded.south.korea$ST[excluded.south.korea$sleeping_time_hour_alternative>6 | excluded.south.korea$sleeping_time_hour_alternative<19] <- "Drop" # Sleeping time > 06:00 & Sleeping time < 19:00
excluded.south.korea$ST <- ifelse(is.na(excluded.south.korea$ST), "Keep", excluded.south.korea$ST) # Sleeping time ≤ 06:00 & Sleeping time ≥ 19:00
table(excluded.south.korea$ST, exclude = NULL)

excluded.south.korea$GT <- NA
excluded.south.korea$GT[excluded.south.korea$getting_up_time_hour_alternative>15 | excluded.south.korea$getting_up_time_hour_alternative<4] <- "Drop" # Getting up time > 15:00 & getting_up time < 04:00
excluded.south.korea$GT <- ifelse(is.na(excluded.south.korea$GT), "Keep", excluded.south.korea$GT) # Getting up time ≤ 15:00 & Getting up time ≥ 04:00
table(excluded.south.korea$GT, exclude = NULL)

excluded.south.korea$SD <- NA
excluded.south.korea$SD[excluded.south.korea$sleep_duration_alternative>18 | excluded.south.korea$sleep_duration_alternative<2] <- "Drop" # "Sleep duration > 18 hrs & Sleep duration < 2 hrs"
excluded.south.korea$SD <- ifelse(is.na(excluded.south.korea$SD), "Keep", excluded.south.korea$SD) # "Sleep duration ≤ 18 hrs & Sleep duration ≥ 2 hrs"
table(excluded.south.korea$SD, exclude = NULL)

excluded.south.korea.table <- xtabs(~ ST + GT + SD, data=excluded.south.korea)
excluded.south.korea.table 
ftable(excluded.south.korea.table)

#### Pooled Data - Grand mean centering ####

issp_2007$PARTY_LR_avg <- mean(issp_2007$PARTY_LR, na.rm = T)
issp_2007$PARTY_LR_sd <- sd(issp_2007$PARTY_LR, na.rm = T)
issp_2007$PARTY_LR_2sd <- issp_2007$PARTY_LR_sd*2
issp_2007$PARTY_LR_beta <- (issp_2007$PARTY_LR-issp_2007$PARTY_LR_avg)/issp_2007$PARTY_LR_2sd

issp_2007$religiosity_avg <- mean(issp_2007$religiosity, na.rm = T)
issp_2007$religiosity_sd <- sd(issp_2007$religiosity, na.rm = T)
issp_2007$religiosity_2sd <- issp_2007$religiosity_sd*2
issp_2007$religiosity_beta <- (issp_2007$religiosity-issp_2007$religiosity_avg)/issp_2007$religiosity_2sd

issp_2007 <- issp_2007 %>%
  group_by(country) %>%
  mutate(group_mean_chronotype = mean(chronotype_interval, na.rm = T),
         group_sd_chronotype = sd(chronotype_interval, na.rm = T))

issp_2007$chronotype_interval_avg <- issp_2007$group_mean_chronotype
issp_2007$chronotype_interval_sd <- issp_2007$group_sd_chronotype
issp_2007$chronotype_interval_2sd <- issp_2007$chronotype_interval_sd*2
issp_2007$chronotype_interval_beta <- (issp_2007$chronotype_interval-issp_2007$chronotype_interval_avg)/issp_2007$chronotype_interval_2sd

issp_2007 <- issp_2007 %>%
  group_by(country) %>%
  mutate(group_mean_sleep_duration = mean(sleep_duration, na.rm = T),
         group_sd_sleep_duration = sd(sleep_duration, na.rm = T))

issp_2007$sleep_duration_avg <- issp_2007$group_mean_sleep_duration
issp_2007$sleep_duration_sd <- issp_2007$group_sd_sleep_duration
issp_2007$sleep_duration_2sd <- issp_2007$sleep_duration_sd*2
issp_2007$sleep_duration_beta <- (issp_2007$sleep_duration-issp_2007$sleep_duration_avg)/issp_2007$sleep_duration_2sd

table(issp_2007$urban_rural_area, exclude = NULL)
issp_2007$urban_rural_area_avg <- mean(issp_2007$urban_rural_area, na.rm = T)
issp_2007$urban_rural_area_sd <- sd(issp_2007$urban_rural_area, na.rm = T)
issp_2007$urban_rural_area_2sd <- issp_2007$urban_rural_area_sd*2
issp_2007$urban_rural_area_beta <- (issp_2007$urban_rural_area-issp_2007$urban_rural_area_avg)/issp_2007$urban_rural_area_2sd

table(issp_2007$age, exclude = NULL)
issp_2007$age_avg <- mean(issp_2007$age, na.rm = T)
issp_2007$age_sd <- sd(issp_2007$age, na.rm = T)
issp_2007$age_2sd <- issp_2007$age_sd*2
issp_2007$age_beta <- (issp_2007$age-issp_2007$age_avg)/issp_2007$age_2sd

table(issp_2007$degree, exclude = NULL)
issp_2007$degree_avg <- mean(issp_2007$degree, na.rm = T)
issp_2007$degree_sd <- sd(issp_2007$degree, na.rm = T)
issp_2007$degree_2sd <- issp_2007$degree_sd*2
issp_2007$degree_beta <- (issp_2007$degree-issp_2007$degree_avg)/issp_2007$degree_2sd

issp_2007$sex_dummy_avg <- mean(issp_2007$sex, na.rm = T)
issp_2007$sex_beta <- (issp_2007$sex-issp_2007$sex_dummy_avg)

issp_2007$political_interest_avg <- mean(issp_2007$political_interest, na.rm = T)
issp_2007$political_interest_sd <- sd(issp_2007$political_interest, na.rm = T)
issp_2007$political_interest_2sd <- issp_2007$political_interest_sd*2
issp_2007$political_interest_beta <- (issp_2007$political_interest-issp_2007$political_interest_avg)/issp_2007$political_interest_2sd

issp_2007$income_avg <- mean(issp_2007$income, na.rm = T)
issp_2007$income_sd <- sd(issp_2007$income, na.rm = T)
issp_2007$income_2sd <- issp_2007$income_sd*2
issp_2007$income_beta <- (issp_2007$income-issp_2007$income_avg)/issp_2007$income_2sd

issp_2007$turnout_dummy_avg <- mean(issp_2007$turnout, na.rm = T)
issp_2007$turnout_beta <- (issp_2007$turnout-issp_2007$turnout_dummy_avg)

issp_2007$polyarchy_avg <- mean(issp_2007$polyarchy, na.rm = T)
issp_2007$polyarchy_sd <- sd(issp_2007$polyarchy, na.rm = T)
issp_2007$polyarchy_2sd <- issp_2007$polyarchy_sd*2
issp_2007$polyarchy_beta <- (issp_2007$polyarchy-issp_2007$polyarchy_avg)/issp_2007$polyarchy_2sd

table(issp_2007$country, issp_2007$mid_field_season, exclude = NULL)
issp_2007$mid_field_season <- ifelse(is.na(issp_2007$mid_field_season) & issp_2007$country=="New Zealand",
                                     "Spring", issp_2007$mid_field_season)
issp_2007$mid_field_season <- ifelse(is.na(issp_2007$mid_field_season) & issp_2007$country=="Greece",
                                     "Winter", issp_2007$mid_field_season)

table(issp_2007$country, exclude = NULL)
issp_2007$election_season <- NA
issp_2007$election_season[issp_2007$country=="Finland" |
                            issp_2007$country=="Ireland" |
                            issp_2007$country=="the Netherlands" |
                            issp_2007$country=="New Zealand" |
                            issp_2007$country=="the Philippines" |
                            issp_2007$country=="Russia" |
                            issp_2007$country=="South Korea"] <- "Spring"
issp_2007$election_season[issp_2007$country=="Greece" |
                            issp_2007$country=="Mexico"] <- "Summer"
table(issp_2007$election_season, exclude = NULL)
issp_2007 <- fastDummies::dummy_cols(issp_2007, select_columns = 'election_season')

table(issp_2007$election_season_Summer, exclude = NULL)
issp_2007$election_season_Summer_avg <- mean(issp_2007$election_season_Summer, na.rm = T)
issp_2007$election_season_Summer_beta <- (issp_2007$election_season_Summer-issp_2007$election_season_Summer_avg)
table(issp_2007$election_season_Summer_beta, exclude = NULL)

table(issp_2007$latitude, exclude = NULL)
issp_2007$latitude_avg <- mean(issp_2007$latitude, na.rm = T)
issp_2007$latitude_sd <- sd(issp_2007$latitude, na.rm = T)
issp_2007$latitude_2sd <- issp_2007$latitude_sd*2
issp_2007$latitude_beta <- (issp_2007$latitude-issp_2007$latitude_avg)/issp_2007$latitude_2sd

table(issp_2007$longitude, exclude = NULL)
issp_2007$longitude_avg <- mean(issp_2007$longitude, na.rm = T)
issp_2007$longitude_sd <- sd(issp_2007$longitude, na.rm = T)
issp_2007$longitude_2sd <- issp_2007$longitude_sd*2
issp_2007$longitude_beta <- (issp_2007$longitude-issp_2007$longitude_avg)/issp_2007$longitude_2sd

table(issp_2007$mean_temperature, exclude = NULL)
issp_2007$mean_temperature_avg <- mean(issp_2007$mean_temperature, na.rm = T)
issp_2007$mean_temperature_sd <- sd(issp_2007$mean_temperature, na.rm = T)
issp_2007$mean_temperature_2sd <- issp_2007$mean_temperature_sd*2
issp_2007$mean_temperature_beta <- (issp_2007$mean_temperature-issp_2007$mean_temperature_avg)/issp_2007$mean_temperature_2sd

table(issp_2007$country, exclude = NULL)
issp_2007$time_duration_election_turnout <- NA
issp_2007$time_duration_election_turnout[issp_2007$country=="Finland"] <- 0
issp_2007$time_duration_election_turnout[issp_2007$country=="Greece"] <- 1
issp_2007$time_duration_election_turnout[issp_2007$country=="Ireland"] <- 1
issp_2007$time_duration_election_turnout[issp_2007$country=="Mexico"] <- 2
issp_2007$time_duration_election_turnout[issp_2007$country=="the Netherlands"] <- 2
issp_2007$time_duration_election_turnout[issp_2007$country=="New Zealand"] <- 2
issp_2007$time_duration_election_turnout[issp_2007$country=="the Philippines"] <- 4
issp_2007$time_duration_election_turnout[issp_2007$country=="Russia"] <- 3
issp_2007$time_duration_election_turnout[issp_2007$country=="South Korea"] <- 1
table(issp_2007$time_duration_election_turnout, exclude = NULL)

issp_2007$time_duration_election_turnout_avg <- mean(issp_2007$time_duration_election_turnout, na.rm = T)
issp_2007$time_duration_election_turnout_sd <- sd(issp_2007$time_duration_election_turnout, na.rm = T)
issp_2007$time_duration_election_turnout_2sd <- issp_2007$time_duration_election_turnout_sd*2
issp_2007$time_duration_election_turnout_beta <- (issp_2007$time_duration_election_turnout-issp_2007$time_duration_election_turnout_avg)/issp_2007$time_duration_election_turnout_2sd
table(issp_2007$time_duration_election_turnout_beta, exclude = NULL)

#### UNREPORTED Descriptive statistics ####

table1::table1(~  turnout + 
                 sleep_duration +  
                 chronotype_interval +  
                 urban_rural_area + 
                 age +  
                 degree +  
                 income +  
                 sex +  
                 political_interest + 
                 PARTY_LR +  
                 religiosity| country, data=issp_2007, digits=2)

#### Main Models ####

issp_2007_turnout_sleep_duration_model1 <- glmer(turnout ~ sleep_duration_beta + chronotype_interval_beta + 
                                                   (1 + sleep_duration_beta + chronotype_interval_beta | country/local_election), 
                                                 data =  issp_2007, family = "binomial", nAGQ=0)
summary(issp_2007_turnout_sleep_duration_model1)
tab_model(issp_2007_turnout_sleep_duration_model1, show.ci = 0.95)

issp_2007_turnout_sleep_duration_model2 <- glmer(turnout ~ sleep_duration_beta + chronotype_interval_beta + 
                                                   urban_rural_area_beta + (1  + sleep_duration_beta + chronotype_interval_beta | country/local_election), 
                                                 data =  issp_2007, family = "binomial", nAGQ=0)
summary(issp_2007_turnout_sleep_duration_model2)
tab_model(issp_2007_turnout_sleep_duration_model2, show.ci = 0.95)

issp_2007_turnout_sleep_duration_model3 <- glmer(turnout ~ sleep_duration_beta + chronotype_interval_beta + 
                                                   urban_rural_area_beta +
                                                   age_beta + degree_beta + income_beta + sex_beta + political_interest_beta + PARTY_LR_beta + 
                                                   (1  + sleep_duration_beta + chronotype_interval_beta | country/local_election), 
                                                 data =  issp_2007, family = "binomial", nAGQ=0)
summary(issp_2007_turnout_sleep_duration_model3)
tab_model(issp_2007_turnout_sleep_duration_model3, show.ci = 0.95)

issp_2007_turnout_sleep_duration_model4 <- glmer(turnout ~ sleep_duration_beta + chronotype_interval_beta + 
                                                   urban_rural_area_beta +
                                                   age_beta + degree_beta + income_beta + sex_beta + political_interest_beta + religiosity_beta + 
                                                   (1  + sleep_duration_beta + chronotype_interval_beta | country/local_election), 
                                                 data =  issp_2007, family = "binomial", nAGQ=0)
summary(issp_2007_turnout_sleep_duration_model4)
tab_model(issp_2007_turnout_sleep_duration_model4, show.ci = 0.95)

####  Online Appendices, Table 14 Multilevel logistic regression models of turnout by sleep duration and chronotype (by including the sample previously excluded). ####

stargazer(issp_2007_turnout_sleep_duration_model1, issp_2007_turnout_sleep_duration_model2, issp_2007_turnout_sleep_duration_model3, issp_2007_turnout_sleep_duration_model4,
          type = "html", title=" ", digits=2, out="Tables/models_issp_2007_turnout_sleep_duration_chronotype_w_excluded.htm",
          model.numbers = F,
          column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),
          covariate.labels = c("Sleep Duration", 
                               "Chronotype",
                               "Urban-Rural Area of Residence", 
                               "Age", 
                               "Education", 
                               "Income", 
                               "Sex: Male (Base: Female)",
                               "Level on Interest in Politics", 
                               "Left-Right Ideological Placement", 
                               "Religious Attendance", 
                               "Intercept"))

####  Manuscript, Table 15 Predicted effects of sleep duration and chronotype on turnout for each country (by including the sample previously excluded). ####

coef(issp_2007_turnout_sleep_duration_model1)
arm::se.coef(issp_2007_turnout_sleep_duration_model1)

coef(issp_2007_turnout_sleep_duration_model1)$country
coef_issp_2007_turnout_sleep_duration_model1 <- data.frame(coef(issp_2007_turnout_sleep_duration_model1)$country)
coef_issp_2007_turnout_sleep_duration_model1 <- subset(coef_issp_2007_turnout_sleep_duration_model1, select = -c(X.Intercept.))
coef_issp_2007_turnout_sleep_duration_model1$country <- row.names(coef_issp_2007_turnout_sleep_duration_model1)
names(coef_issp_2007_turnout_sleep_duration_model1)[names(coef_issp_2007_turnout_sleep_duration_model1)=="sleep_duration_beta"] <- "coef_sleep_duration_beta"
names(coef_issp_2007_turnout_sleep_duration_model1)[names(coef_issp_2007_turnout_sleep_duration_model1)=="chronotype_interval_beta"] <- "coef_chronotype_interval_beta"


se_coef_issp_2007_turnout_sleep_duration_model1 <- data.frame(arm::se.coef(issp_2007_turnout_sleep_duration_model1)$country)
se_coef_issp_2007_turnout_sleep_duration_model1 <- subset(se_coef_issp_2007_turnout_sleep_duration_model1, select = -c(X.Intercept.))
se_coef_issp_2007_turnout_sleep_duration_model1$country <- row.names(se_coef_issp_2007_turnout_sleep_duration_model1)
names(se_coef_issp_2007_turnout_sleep_duration_model1)[names(se_coef_issp_2007_turnout_sleep_duration_model1)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_coef_issp_2007_turnout_sleep_duration_model1)[names(se_coef_issp_2007_turnout_sleep_duration_model1)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

coef_se_issp_2007_turnout_sleep_duration_model1 <- merge(coef_issp_2007_turnout_sleep_duration_model1, se_coef_issp_2007_turnout_sleep_duration_model1, by="country")

coef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_sleep_duration_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_sleep_duration_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_sleep_duration_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_sleep_duration_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_sleep_duration_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_sleep_duration_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta

coef_se_issp_2007_turnout_sleep_duration_model1$chro.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_chronotype_interval_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model1$chro.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_chronotype_interval_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model1$chro.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_chronotype_interval_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model1$chro.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_chronotype_interval_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model1$chro.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_chronotype_interval_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model1$chro.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model1$coef_chronotype_interval_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta

View(coef_se_issp_2007_turnout_sleep_duration_model1)

coef(issp_2007_turnout_sleep_duration_model2)
arm::se.coef(issp_2007_turnout_sleep_duration_model2)

coef(issp_2007_turnout_sleep_duration_model2)$country
coef_issp_2007_turnout_sleep_duration_model2 <- data.frame(coef(issp_2007_turnout_sleep_duration_model2)$country)
coef_issp_2007_turnout_sleep_duration_model2 <- subset(coef_issp_2007_turnout_sleep_duration_model2, select = c("sleep_duration_beta", "chronotype_interval_beta"))
coef_issp_2007_turnout_sleep_duration_model2$country <- row.names(coef_issp_2007_turnout_sleep_duration_model2)
names(coef_issp_2007_turnout_sleep_duration_model2)[names(coef_issp_2007_turnout_sleep_duration_model2)=="sleep_duration_beta"] <- "coef_sleep_duration_beta"
names(coef_issp_2007_turnout_sleep_duration_model2)[names(coef_issp_2007_turnout_sleep_duration_model2)=="chronotype_interval_beta"] <- "coef_chronotype_interval_beta"


se_coef_issp_2007_turnout_sleep_duration_model2 <- data.frame(arm::se.coef(issp_2007_turnout_sleep_duration_model2)$country)
se_coef_issp_2007_turnout_sleep_duration_model2 <- subset(se_coef_issp_2007_turnout_sleep_duration_model2, select = c("sleep_duration_beta", "chronotype_interval_beta"))
se_coef_issp_2007_turnout_sleep_duration_model2$country <- row.names(se_coef_issp_2007_turnout_sleep_duration_model2)
names(se_coef_issp_2007_turnout_sleep_duration_model2)[names(se_coef_issp_2007_turnout_sleep_duration_model2)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_coef_issp_2007_turnout_sleep_duration_model2)[names(se_coef_issp_2007_turnout_sleep_duration_model2)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

coef_se_issp_2007_turnout_sleep_duration_model2 <- merge(coef_issp_2007_turnout_sleep_duration_model2, se_coef_issp_2007_turnout_sleep_duration_model2, by="country")

coef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_sleep_duration_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_sleep_duration_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_sleep_duration_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_sleep_duration_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_sleep_duration_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_sleep_duration_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta

coef_se_issp_2007_turnout_sleep_duration_model2$chro.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_chronotype_interval_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model2$chro.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_chronotype_interval_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model2$chro.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_chronotype_interval_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model2$chro.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_chronotype_interval_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model2$chro.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_chronotype_interval_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model2$chro.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model2$coef_chronotype_interval_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta

View(coef_se_issp_2007_turnout_sleep_duration_model2)

coef(issp_2007_turnout_sleep_duration_model3)
arm::se.coef(issp_2007_turnout_sleep_duration_model3)

coef(issp_2007_turnout_sleep_duration_model3)$country
coef_issp_2007_turnout_sleep_duration_model3 <- data.frame(coef(issp_2007_turnout_sleep_duration_model3)$country)
coef_issp_2007_turnout_sleep_duration_model3 <- subset(coef_issp_2007_turnout_sleep_duration_model3, select = c("sleep_duration_beta", "chronotype_interval_beta"))
coef_issp_2007_turnout_sleep_duration_model3$country <- row.names(coef_issp_2007_turnout_sleep_duration_model3)
names(coef_issp_2007_turnout_sleep_duration_model3)[names(coef_issp_2007_turnout_sleep_duration_model3)=="sleep_duration_beta"] <- "coef_sleep_duration_beta"
names(coef_issp_2007_turnout_sleep_duration_model3)[names(coef_issp_2007_turnout_sleep_duration_model3)=="chronotype_interval_beta"] <- "coef_chronotype_interval_beta"


se_coef_issp_2007_turnout_sleep_duration_model3 <- data.frame(arm::se.coef(issp_2007_turnout_sleep_duration_model3)$country)
se_coef_issp_2007_turnout_sleep_duration_model3 <- subset(se_coef_issp_2007_turnout_sleep_duration_model3, select = c("sleep_duration_beta", "chronotype_interval_beta"))
se_coef_issp_2007_turnout_sleep_duration_model3$country <- row.names(se_coef_issp_2007_turnout_sleep_duration_model3)
names(se_coef_issp_2007_turnout_sleep_duration_model3)[names(se_coef_issp_2007_turnout_sleep_duration_model3)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_coef_issp_2007_turnout_sleep_duration_model3)[names(se_coef_issp_2007_turnout_sleep_duration_model3)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

coef_se_issp_2007_turnout_sleep_duration_model3 <- merge(coef_issp_2007_turnout_sleep_duration_model3, se_coef_issp_2007_turnout_sleep_duration_model3, by="country")

coef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_sleep_duration_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_sleep_duration_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_sleep_duration_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_sleep_duration_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_sleep_duration_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_sleep_duration_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta

coef_se_issp_2007_turnout_sleep_duration_model3$chro.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_chronotype_interval_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model3$chro.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_chronotype_interval_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model3$chro.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_chronotype_interval_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model3$chro.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_chronotype_interval_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model3$chro.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_chronotype_interval_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model3$chro.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model3$coef_chronotype_interval_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta

View(coef_se_issp_2007_turnout_sleep_duration_model3)

coef(issp_2007_turnout_sleep_duration_model4)
arm::se.coef(issp_2007_turnout_sleep_duration_model4)

coef(issp_2007_turnout_sleep_duration_model4)$country
coef_issp_2007_turnout_sleep_duration_model4 <- data.frame(coef(issp_2007_turnout_sleep_duration_model4)$country)
coef_issp_2007_turnout_sleep_duration_model4 <- subset(coef_issp_2007_turnout_sleep_duration_model4, select = c("sleep_duration_beta", "chronotype_interval_beta"))
coef_issp_2007_turnout_sleep_duration_model4$country <- row.names(coef_issp_2007_turnout_sleep_duration_model4)
names(coef_issp_2007_turnout_sleep_duration_model4)[names(coef_issp_2007_turnout_sleep_duration_model4)=="sleep_duration_beta"] <- "coef_sleep_duration_beta"
names(coef_issp_2007_turnout_sleep_duration_model4)[names(coef_issp_2007_turnout_sleep_duration_model4)=="chronotype_interval_beta"] <- "coef_chronotype_interval_beta"


se_coef_issp_2007_turnout_sleep_duration_model4 <- data.frame(arm::se.coef(issp_2007_turnout_sleep_duration_model4)$country)
se_coef_issp_2007_turnout_sleep_duration_model4 <- subset(se_coef_issp_2007_turnout_sleep_duration_model4, select = c("sleep_duration_beta", "chronotype_interval_beta"))
se_coef_issp_2007_turnout_sleep_duration_model4$country <- row.names(se_coef_issp_2007_turnout_sleep_duration_model4)
names(se_coef_issp_2007_turnout_sleep_duration_model4)[names(se_coef_issp_2007_turnout_sleep_duration_model4)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_coef_issp_2007_turnout_sleep_duration_model4)[names(se_coef_issp_2007_turnout_sleep_duration_model4)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

coef_se_issp_2007_turnout_sleep_duration_model4 <- merge(coef_issp_2007_turnout_sleep_duration_model4, se_coef_issp_2007_turnout_sleep_duration_model4, by="country")

coef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_sleep_duration_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_sleep_duration_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_sleep_duration_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_sleep_duration_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_sleep_duration_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
coef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_sleep_duration_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta

coef_se_issp_2007_turnout_sleep_duration_model4$chro.cl.90 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_chronotype_interval_beta - 1.645*coef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model4$chro.cu.90 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_chronotype_interval_beta + 1.645*coef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model4$chro.cl.95 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_chronotype_interval_beta - 1.96*coef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model4$chro.cu.95 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_chronotype_interval_beta + 1.96*coef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model4$chro.cl.99 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_chronotype_interval_beta - 2.576*coef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
coef_se_issp_2007_turnout_sleep_duration_model4$chro.cu.99 <- coef_se_issp_2007_turnout_sleep_duration_model4$coef_chronotype_interval_beta + 2.576*coef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta

View(coef_se_issp_2007_turnout_sleep_duration_model4)

####  Online Appendices, Table 16 Random effects of sleep duration and chronotype on turnout for each country (by including the sample previously excluded). ####

ranef(issp_2007_turnout_sleep_duration_model1)
arm::se.ranef(issp_2007_turnout_sleep_duration_model1)

ranef(issp_2007_turnout_sleep_duration_model1)$country
ranef_issp_2007_turnout_sleep_duration_model1 <- data.frame(ranef(issp_2007_turnout_sleep_duration_model1)$country)
ranef_issp_2007_turnout_sleep_duration_model1 <- subset(ranef_issp_2007_turnout_sleep_duration_model1, select = -c(X.Intercept.))
ranef_issp_2007_turnout_sleep_duration_model1$country <- row.names(ranef_issp_2007_turnout_sleep_duration_model1)
names(ranef_issp_2007_turnout_sleep_duration_model1)[names(ranef_issp_2007_turnout_sleep_duration_model1)=="sleep_duration_beta"] <- "ranef_sleep_duration_beta"
names(ranef_issp_2007_turnout_sleep_duration_model1)[names(ranef_issp_2007_turnout_sleep_duration_model1)=="chronotype_interval_beta"] <- "ranef_chronotype_interval_beta"


se_ranef_issp_2007_turnout_sleep_duration_model1 <- data.frame(arm::se.ranef(issp_2007_turnout_sleep_duration_model1)$country)
se_ranef_issp_2007_turnout_sleep_duration_model1 <- subset(se_ranef_issp_2007_turnout_sleep_duration_model1, select = -c(X.Intercept.))
se_ranef_issp_2007_turnout_sleep_duration_model1$country <- row.names(se_ranef_issp_2007_turnout_sleep_duration_model1)
names(se_ranef_issp_2007_turnout_sleep_duration_model1)[names(se_ranef_issp_2007_turnout_sleep_duration_model1)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_ranef_issp_2007_turnout_sleep_duration_model1)[names(se_ranef_issp_2007_turnout_sleep_duration_model1)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

ranef_se_issp_2007_turnout_sleep_duration_model1 <- merge(ranef_issp_2007_turnout_sleep_duration_model1, se_ranef_issp_2007_turnout_sleep_duration_model1, by="country")

ranef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_sleep_duration_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_sleep_duration_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_sleep_duration_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_sleep_duration_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_sleep_duration_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$slp.dur.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_sleep_duration_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model1$se_sleep_duration_beta

ranef_se_issp_2007_turnout_sleep_duration_model1$chro.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_chronotype_interval_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$chro.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_chronotype_interval_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$chro.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_chronotype_interval_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$chro.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_chronotype_interval_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$chro.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_chronotype_interval_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model1$chro.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model1$ranef_chronotype_interval_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model1$se_chronotype_interval_beta

View(ranef_se_issp_2007_turnout_sleep_duration_model1)

ranef(issp_2007_turnout_sleep_duration_model2)
arm::se.ranef(issp_2007_turnout_sleep_duration_model2)

ranef(issp_2007_turnout_sleep_duration_model2)$country
ranef_issp_2007_turnout_sleep_duration_model2 <- data.frame(ranef(issp_2007_turnout_sleep_duration_model2)$country)
ranef_issp_2007_turnout_sleep_duration_model2 <- subset(ranef_issp_2007_turnout_sleep_duration_model2, select = c("sleep_duration_beta", "chronotype_interval_beta"))
ranef_issp_2007_turnout_sleep_duration_model2$country <- row.names(ranef_issp_2007_turnout_sleep_duration_model2)
names(ranef_issp_2007_turnout_sleep_duration_model2)[names(ranef_issp_2007_turnout_sleep_duration_model2)=="sleep_duration_beta"] <- "ranef_sleep_duration_beta"
names(ranef_issp_2007_turnout_sleep_duration_model2)[names(ranef_issp_2007_turnout_sleep_duration_model2)=="chronotype_interval_beta"] <- "ranef_chronotype_interval_beta"


se_ranef_issp_2007_turnout_sleep_duration_model2 <- data.frame(arm::se.ranef(issp_2007_turnout_sleep_duration_model2)$country)
se_ranef_issp_2007_turnout_sleep_duration_model2 <- subset(se_ranef_issp_2007_turnout_sleep_duration_model2, select = c("sleep_duration_beta", "chronotype_interval_beta"))
se_ranef_issp_2007_turnout_sleep_duration_model2$country <- row.names(se_ranef_issp_2007_turnout_sleep_duration_model2)
names(se_ranef_issp_2007_turnout_sleep_duration_model2)[names(se_ranef_issp_2007_turnout_sleep_duration_model2)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_ranef_issp_2007_turnout_sleep_duration_model2)[names(se_ranef_issp_2007_turnout_sleep_duration_model2)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

ranef_se_issp_2007_turnout_sleep_duration_model2 <- merge(ranef_issp_2007_turnout_sleep_duration_model2, se_ranef_issp_2007_turnout_sleep_duration_model2, by="country")

ranef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_sleep_duration_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_sleep_duration_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_sleep_duration_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_sleep_duration_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_sleep_duration_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$slp.dur.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_sleep_duration_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model2$se_sleep_duration_beta

ranef_se_issp_2007_turnout_sleep_duration_model2$chro.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_chronotype_interval_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$chro.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_chronotype_interval_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$chro.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_chronotype_interval_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$chro.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_chronotype_interval_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$chro.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_chronotype_interval_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model2$chro.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model2$ranef_chronotype_interval_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model2$se_chronotype_interval_beta

View(ranef_se_issp_2007_turnout_sleep_duration_model2)

ranef(issp_2007_turnout_sleep_duration_model3)
arm::se.ranef(issp_2007_turnout_sleep_duration_model3)

ranef(issp_2007_turnout_sleep_duration_model3)$country
ranef_issp_2007_turnout_sleep_duration_model3 <- data.frame(ranef(issp_2007_turnout_sleep_duration_model3)$country)
ranef_issp_2007_turnout_sleep_duration_model3 <- subset(ranef_issp_2007_turnout_sleep_duration_model3, select = c("sleep_duration_beta", "chronotype_interval_beta"))
ranef_issp_2007_turnout_sleep_duration_model3$country <- row.names(ranef_issp_2007_turnout_sleep_duration_model3)
names(ranef_issp_2007_turnout_sleep_duration_model3)[names(ranef_issp_2007_turnout_sleep_duration_model3)=="sleep_duration_beta"] <- "ranef_sleep_duration_beta"
names(ranef_issp_2007_turnout_sleep_duration_model3)[names(ranef_issp_2007_turnout_sleep_duration_model3)=="chronotype_interval_beta"] <- "ranef_chronotype_interval_beta"


se_ranef_issp_2007_turnout_sleep_duration_model3 <- data.frame(arm::se.ranef(issp_2007_turnout_sleep_duration_model3)$country)
se_ranef_issp_2007_turnout_sleep_duration_model3 <- subset(se_ranef_issp_2007_turnout_sleep_duration_model3, select = c("sleep_duration_beta", "chronotype_interval_beta"))
se_ranef_issp_2007_turnout_sleep_duration_model3$country <- row.names(se_ranef_issp_2007_turnout_sleep_duration_model3)
names(se_ranef_issp_2007_turnout_sleep_duration_model3)[names(se_ranef_issp_2007_turnout_sleep_duration_model3)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_ranef_issp_2007_turnout_sleep_duration_model3)[names(se_ranef_issp_2007_turnout_sleep_duration_model3)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

ranef_se_issp_2007_turnout_sleep_duration_model3 <- merge(ranef_issp_2007_turnout_sleep_duration_model3, se_ranef_issp_2007_turnout_sleep_duration_model3, by="country")

ranef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_sleep_duration_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_sleep_duration_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_sleep_duration_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_sleep_duration_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_sleep_duration_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$slp.dur.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_sleep_duration_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model3$se_sleep_duration_beta

ranef_se_issp_2007_turnout_sleep_duration_model3$chro.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_chronotype_interval_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$chro.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_chronotype_interval_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$chro.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_chronotype_interval_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$chro.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_chronotype_interval_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$chro.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_chronotype_interval_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model3$chro.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model3$ranef_chronotype_interval_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model3$se_chronotype_interval_beta

View(ranef_se_issp_2007_turnout_sleep_duration_model3)

ranef(issp_2007_turnout_sleep_duration_model4)
arm::se.ranef(issp_2007_turnout_sleep_duration_model4)

ranef(issp_2007_turnout_sleep_duration_model4)$country
ranef_issp_2007_turnout_sleep_duration_model4 <- data.frame(ranef(issp_2007_turnout_sleep_duration_model4)$country)
ranef_issp_2007_turnout_sleep_duration_model4 <- subset(ranef_issp_2007_turnout_sleep_duration_model4, select = c("sleep_duration_beta", "chronotype_interval_beta"))
ranef_issp_2007_turnout_sleep_duration_model4$country <- row.names(ranef_issp_2007_turnout_sleep_duration_model4)
names(ranef_issp_2007_turnout_sleep_duration_model4)[names(ranef_issp_2007_turnout_sleep_duration_model4)=="sleep_duration_beta"] <- "ranef_sleep_duration_beta"
names(ranef_issp_2007_turnout_sleep_duration_model4)[names(ranef_issp_2007_turnout_sleep_duration_model4)=="chronotype_interval_beta"] <- "ranef_chronotype_interval_beta"


se_ranef_issp_2007_turnout_sleep_duration_model4 <- data.frame(arm::se.ranef(issp_2007_turnout_sleep_duration_model4)$country)
se_ranef_issp_2007_turnout_sleep_duration_model4 <- subset(se_ranef_issp_2007_turnout_sleep_duration_model4, select = c("sleep_duration_beta", "chronotype_interval_beta"))
se_ranef_issp_2007_turnout_sleep_duration_model4$country <- row.names(se_ranef_issp_2007_turnout_sleep_duration_model4)
names(se_ranef_issp_2007_turnout_sleep_duration_model4)[names(se_ranef_issp_2007_turnout_sleep_duration_model4)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_ranef_issp_2007_turnout_sleep_duration_model4)[names(se_ranef_issp_2007_turnout_sleep_duration_model4)=="chronotype_interval_beta"] <- "se_chronotype_interval_beta"

ranef_se_issp_2007_turnout_sleep_duration_model4 <- merge(ranef_issp_2007_turnout_sleep_duration_model4, se_ranef_issp_2007_turnout_sleep_duration_model4, by="country")

ranef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_sleep_duration_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_sleep_duration_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_sleep_duration_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_sleep_duration_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_sleep_duration_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$slp.dur.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_sleep_duration_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model4$se_sleep_duration_beta

ranef_se_issp_2007_turnout_sleep_duration_model4$chro.cl.90 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_chronotype_interval_beta - 1.645*ranef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$chro.cu.90 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_chronotype_interval_beta + 1.645*ranef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$chro.cl.95 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_chronotype_interval_beta - 1.96*ranef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$chro.cu.95 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_chronotype_interval_beta + 1.96*ranef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$chro.cl.99 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_chronotype_interval_beta - 2.576*ranef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta
ranef_se_issp_2007_turnout_sleep_duration_model4$chro.cu.99 <- ranef_se_issp_2007_turnout_sleep_duration_model4$ranef_chronotype_interval_beta + 2.576*ranef_se_issp_2007_turnout_sleep_duration_model4$se_chronotype_interval_beta

View(ranef_se_issp_2007_turnout_sleep_duration_model4)

#### Online Appendices, Table 17 Multilevel logistic regression models of turnout by sleep duration and sleep duration squared with controls (by including the sample previously excluded). ####

issp_2007$sleep_duration_power2 <- issp_2007$sleep_duration^2
hist(issp_2007$sleep_duration_power2)

issp_2007 <- issp_2007 %>%
  group_by(country) %>%
  mutate(group_mean_sleep_duration_power2 = mean(sleep_duration_power2, na.rm = T),
         group_sd_sleep_duration_power2 = sd(sleep_duration_power2, na.rm = T))

issp_2007$sleep_duration_power2_avg <- issp_2007$group_mean_sleep_duration_power2
issp_2007$sleep_duration_power2_sd <- issp_2007$group_sd_sleep_duration_power2
issp_2007$sleep_duration_power2_2sd <- issp_2007$sleep_duration_power2_sd*2
issp_2007$sleep_duration_power2_beta <- (issp_2007$sleep_duration_power2-issp_2007$sleep_duration_power2_avg)/issp_2007$sleep_duration_power2_2sd
hist(issp_2007$sleep_duration_power2_beta)

poly.issp.m1 <- glmer(turnout ~ sleep_duration_beta + sleep_duration_power2_beta + chronotype_interval_beta +
                        urban_rural_area_beta +
                        age_beta + degree_beta + income_beta + sex_beta + political_interest_beta + religiosity_beta +
                        (1  + sleep_duration_beta + sleep_duration_power2_beta + chronotype_interval_beta | country/local_election) +
                        (1 + sleep_duration_beta + sleep_duration_power2_beta + chronotype_interval_beta | mid_field_season/day_off), 
                      data =  issp_2007, family = "binomial", nAGQ=0,
                      control=glmerControl(optCtrl=list(maxfun=1e5)))
summary(poly.issp.m1)
tab_model(poly.issp.m1, show.ci = 0.95)

stargazer(poly.issp.m1,
          type = "html", title=" ", digits=2, out="Tables/models_polynomial_w_excluded.htm",
          model.numbers = F,
          column.labels = c("Model 1"),
          covariate.labels = c("Sleep Duration",
                               "Sleep Duration Squared",
                               "Chronotype",
                               "Urban-Rural Area of Residence", 
                               "Age", 
                               "Education", 
                               "Income", 
                               "Sex: Male (Base: Female)",
                               "Level on Interest in Politics", 
                               "Religious Attendance", 
                               "Intercept"))

#### Online Appendices, Table 18 Predicted effects of sleep duration and sleep duration squared on turnout for each country (by including the sample previously excluded). ####

coef(poly.issp.m1)
arm::se.coef(poly.issp.m1)

coef(poly.issp.m1)$country
coef_poly.issp.m1 <- data.frame(coef(poly.issp.m1)$country)
coef_poly.issp.m1 <- subset(coef_poly.issp.m1, select = c("sleep_duration_beta", "sleep_duration_power2_beta"))
coef_poly.issp.m1$country <- row.names(coef_poly.issp.m1)
names(coef_poly.issp.m1)[names(coef_poly.issp.m1)=="sleep_duration_beta"] <- "coef_sleep_duration_beta"
names(coef_poly.issp.m1)[names(coef_poly.issp.m1)=="sleep_duration_power2_beta"] <- "coef_sleep_duration_power2_beta"


se_coef_poly.issp.m1 <- data.frame(arm::se.coef(poly.issp.m1)$country)
se_coef_poly.issp.m1 <- subset(se_coef_poly.issp.m1, select = c("sleep_duration_beta", "sleep_duration_power2_beta"))
se_coef_poly.issp.m1$country <- row.names(se_coef_poly.issp.m1)
names(se_coef_poly.issp.m1)[names(se_coef_poly.issp.m1)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_coef_poly.issp.m1)[names(se_coef_poly.issp.m1)=="sleep_duration_power2_beta"] <- "se_sleep_duration_power2_beta"

coef_se_poly.issp.m1 <- merge(coef_poly.issp.m1, se_coef_poly.issp.m1, by="country")

coef_se_poly.issp.m1$slp.dur.cl.90 <- coef_se_poly.issp.m1$coef_sleep_duration_beta - 1.645*coef_se_poly.issp.m1$se_sleep_duration_beta
coef_se_poly.issp.m1$slp.dur.cu.90 <- coef_se_poly.issp.m1$coef_sleep_duration_beta + 1.645*coef_se_poly.issp.m1$se_sleep_duration_beta
coef_se_poly.issp.m1$slp.dur.cl.95 <- coef_se_poly.issp.m1$coef_sleep_duration_beta - 1.96*coef_se_poly.issp.m1$se_sleep_duration_beta
coef_se_poly.issp.m1$slp.dur.cu.95 <- coef_se_poly.issp.m1$coef_sleep_duration_beta + 1.96*coef_se_poly.issp.m1$se_sleep_duration_beta
coef_se_poly.issp.m1$slp.dur.cl.99 <- coef_se_poly.issp.m1$coef_sleep_duration_beta - 2.576*coef_se_poly.issp.m1$se_sleep_duration_beta
coef_se_poly.issp.m1$slp.dur.cu.99 <- coef_se_poly.issp.m1$coef_sleep_duration_beta + 2.576*coef_se_poly.issp.m1$se_sleep_duration_beta

coef_se_poly.issp.m1$slp.dur2.cl.90 <- coef_se_poly.issp.m1$coef_sleep_duration_power2_beta - 1.645*coef_se_poly.issp.m1$se_sleep_duration_power2_beta
coef_se_poly.issp.m1$slp.dur2.cu.90 <- coef_se_poly.issp.m1$coef_sleep_duration_power2_beta + 1.645*coef_se_poly.issp.m1$se_sleep_duration_power2_beta
coef_se_poly.issp.m1$slp.dur2.cl.95 <- coef_se_poly.issp.m1$coef_sleep_duration_power2_beta - 1.96*coef_se_poly.issp.m1$se_sleep_duration_power2_beta
coef_se_poly.issp.m1$slp.dur2.cu.95 <- coef_se_poly.issp.m1$coef_sleep_duration_power2_beta + 1.96*coef_se_poly.issp.m1$se_sleep_duration_power2_beta
coef_se_poly.issp.m1$slp.dur2.cl.99 <- coef_se_poly.issp.m1$coef_sleep_duration_power2_beta - 2.576*coef_se_poly.issp.m1$se_sleep_duration_power2_beta
coef_se_poly.issp.m1$slp.dur2.cu.99 <- coef_se_poly.issp.m1$coef_sleep_duration_power2_beta + 2.576*coef_se_poly.issp.m1$se_sleep_duration_power2_beta

View(coef_se_poly.issp.m1)

#### Online Appendices, Table 19 Random effects of sleep duration and sleep duration squared on turnout for each country (by including the sample previously excluded). ####

ranef(poly.issp.m1)
arm::se.ranef(poly.issp.m1)

ranef(poly.issp.m1)$country
ranef_poly.issp.m1 <- data.frame(ranef(poly.issp.m1)$country)
ranef_poly.issp.m1 <- subset(ranef_poly.issp.m1, select = c("sleep_duration_beta", "sleep_duration_power2_beta"))
ranef_poly.issp.m1$country <- row.names(ranef_poly.issp.m1)
names(ranef_poly.issp.m1)[names(ranef_poly.issp.m1)=="sleep_duration_beta"] <- "ranef_sleep_duration_beta"
names(ranef_poly.issp.m1)[names(ranef_poly.issp.m1)=="sleep_duration_power2_beta"] <- "ranef_sleep_duration_power2_beta"


se_ranef_poly.issp.m1 <- data.frame(arm::se.ranef(poly.issp.m1)$country)
se_ranef_poly.issp.m1 <- subset(se_ranef_poly.issp.m1, select = c("sleep_duration_beta", "sleep_duration_power2_beta"))
se_ranef_poly.issp.m1$country <- row.names(se_ranef_poly.issp.m1)
names(se_ranef_poly.issp.m1)[names(se_ranef_poly.issp.m1)=="sleep_duration_beta"] <- "se_sleep_duration_beta"
names(se_ranef_poly.issp.m1)[names(se_ranef_poly.issp.m1)=="sleep_duration_power2_beta"] <- "se_sleep_duration_power2_beta"

ranef_se_poly.issp.m1 <- merge(ranef_poly.issp.m1, se_ranef_poly.issp.m1, by="country")

ranef_se_poly.issp.m1$slp.dur.cl.90 <- ranef_se_poly.issp.m1$ranef_sleep_duration_beta - 1.645*ranef_se_poly.issp.m1$se_sleep_duration_beta
ranef_se_poly.issp.m1$slp.dur.cu.90 <- ranef_se_poly.issp.m1$ranef_sleep_duration_beta + 1.645*ranef_se_poly.issp.m1$se_sleep_duration_beta
ranef_se_poly.issp.m1$slp.dur.cl.95 <- ranef_se_poly.issp.m1$ranef_sleep_duration_beta - 1.96*ranef_se_poly.issp.m1$se_sleep_duration_beta
ranef_se_poly.issp.m1$slp.dur.cu.95 <- ranef_se_poly.issp.m1$ranef_sleep_duration_beta + 1.96*ranef_se_poly.issp.m1$se_sleep_duration_beta
ranef_se_poly.issp.m1$slp.dur.cl.99 <- ranef_se_poly.issp.m1$ranef_sleep_duration_beta - 2.576*ranef_se_poly.issp.m1$se_sleep_duration_beta
ranef_se_poly.issp.m1$slp.dur.cu.99 <- ranef_se_poly.issp.m1$ranef_sleep_duration_beta + 2.576*ranef_se_poly.issp.m1$se_sleep_duration_beta

ranef_se_poly.issp.m1$slp.dur2.cl.90 <- ranef_se_poly.issp.m1$ranef_sleep_duration_power2_beta - 1.645*ranef_se_poly.issp.m1$se_sleep_duration_power2_beta
ranef_se_poly.issp.m1$slp.dur2.cu.90 <- ranef_se_poly.issp.m1$ranef_sleep_duration_power2_beta + 1.645*ranef_se_poly.issp.m1$se_sleep_duration_power2_beta
ranef_se_poly.issp.m1$slp.dur2.cl.95 <- ranef_se_poly.issp.m1$ranef_sleep_duration_power2_beta - 1.96*ranef_se_poly.issp.m1$se_sleep_duration_power2_beta
ranef_se_poly.issp.m1$slp.dur2.cu.95 <- ranef_se_poly.issp.m1$ranef_sleep_duration_power2_beta + 1.96*ranef_se_poly.issp.m1$se_sleep_duration_power2_beta
ranef_se_poly.issp.m1$slp.dur2.cl.99 <- ranef_se_poly.issp.m1$ranef_sleep_duration_power2_beta - 2.576*ranef_se_poly.issp.m1$se_sleep_duration_power2_beta
ranef_se_poly.issp.m1$slp.dur2.cu.99 <- ranef_se_poly.issp.m1$ranef_sleep_duration_power2_beta + 2.576*ranef_se_poly.issp.m1$se_sleep_duration_power2_beta

View(ranef_se_poly.issp.m1)

#### Online Appendices, Fig. 2. Fitted values of turnout by sleep duration and sleep duration squared by including the sample previously excluded). ####

number_ticks <- function(n) {function(limits) pretty(limits, n)}

table(issp_2007$country, exclude = NULL)

issp_2007$predicted_finland <- predict(poly.issp.m1, 
                                       newdata=transform(issp_2007, 
                                                         urban_rural_area_beta = 0,
                                                         age_beta  = 0, 
                                                         degree_beta  = 0,
                                                         income_beta   = 0,
                                                         sex_beta   = 0,
                                                         political_interest_beta   = 0,
                                                         religiosity_beta = 0,
                                                         country="Finland",
                                                         local_election=0,
                                                         mid_field_season=0,
                                                         day_off=0), type = "response", allow.new.levels=T)

finland <- subset(issp_2007, subset = country == "Finland", select = c("sleep_duration", "predicted_finland"))

describe(finland$predicted_finland)
table(finland$predicted_finland)
describe(finland$sleep_duration)

finland_plot <- ggplot(finland, aes(x=sleep_duration, y=predicted_finland)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(0, 2.91, 5.76, 8.61, 11.46, 14.31, 17.16, 20.01, 22.75),
                     labels = c("0\n(Min.)", "2.91\n(μ-2SD)", "5.76\n(μ-1SD)", "8.61\n(μ)", "11.46\n(μ+1SD)", "14.31\n(μ+2SD)", "17.16\n(μ+3SD)", "20.01\n(μ+4SD)", "22.75\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.70, 0.95), breaks = seq(0.70,0.95, by = .05)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
finland_plot

finland_plot <- ggExtra::ggMarginal(finland_plot, type = "boxplot")

finland_plot

issp_2007$predicted_greece <- predict(poly.issp.m1, 
                                      newdata=transform(issp_2007, 
                                                        urban_rural_area_beta = 0,
                                                        age_beta  = 0, 
                                                        degree_beta  = 0,
                                                        income_beta   = 0,
                                                        sex_beta   = 0,
                                                        political_interest_beta   = 0,
                                                        religiosity_beta = 0,
                                                        country="Greece",
                                                        local_election=0,
                                                        mid_field_season=0,
                                                        day_off=0), type = "response", allow.new.levels=T)

greece <- subset(issp_2007, subset = country == "Greece", select = c("sleep_duration", "predicted_greece"))

describe(greece$predicted_greece)
table(greece$predicted_greece)
describe(greece$sleep_duration)

greece_plot <- ggplot(greece, aes(x=sleep_duration, y=predicted_greece)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(0, 2.23, 5.1, 7.9, 10.7, 13.5, 16.3, 20),
                     labels = c("0\n(Min.)", "2.26\n(μ-4SD)", "5.1\n(μ-2SD)",  "7.9\n(μ)",  "10.7\n(μ+2SD)",  "13.5\n(μ+4SD)", "16.3\n(μ+6SD)", "20\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.87, 0.92), breaks = seq(0.87, 0.92, by = .01)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
greece_plot

greece_plot <- ggExtra::ggMarginal(greece_plot, type = "boxplot")

greece_plot

issp_2007$predicted_ireland <- predict(poly.issp.m1, 
                                       newdata=transform(issp_2007, 
                                                         urban_rural_area_beta = 0,
                                                         age_beta  = 0, 
                                                         degree_beta  = 0,
                                                         income_beta   = 0,
                                                         sex_beta   = 0,
                                                         political_interest_beta   = 0,
                                                         religiosity_beta = 0,
                                                         country="Ireland",
                                                         local_election=0,
                                                         mid_field_season=0,
                                                         day_off=0), type = "response", allow.new.levels=T)

ireland <- subset(issp_2007, subset = country == "Ireland", select = c("sleep_duration", "predicted_ireland"))

describe(ireland$predicted_ireland)
table(ireland$predicted_ireland)
describe(ireland$sleep_duration)

ireland_plot <- ggplot(ireland, aes(x=sleep_duration, y=predicted_ireland)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(2, 3.78, 6.25, 8.72, 11.19, 13.66, 16.13, 18.6, 21.07, 22.92),
                     labels = c("2\n(Min.)", "3.78\n(μ-2SD)", "6.25\n(μ-1SD)", "8.72\n(μ)", "11.19\n(μ+1SD)", "13.66\n(μ+2SD)", "16.13\n(μ+3SD)", "18.6\n(μ+4SD)", "21.07\n(μ+5SD)", "22.92\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.45, 0.90), breaks = seq(0.45, 0.90, by = .05)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
ireland_plot

ireland_plot <- ggExtra::ggMarginal(ireland_plot, type = "boxplot")

ireland_plot

issp_2007$predicted_mexico <- predict(poly.issp.m1, 
                                      newdata=transform(issp_2007, 
                                                        urban_rural_area_beta = 0,
                                                        age_beta  = 0, 
                                                        degree_beta  = 0,
                                                        income_beta   = 0,
                                                        sex_beta   = 0,
                                                        political_interest_beta   = 0,
                                                        religiosity_beta = 0,
                                                        country="Mexico",
                                                        local_election=0,
                                                        mid_field_season=0,
                                                        day_off=0), type = "response", allow.new.levels=T)

mexico <- subset(issp_2007, subset = country == "Mexico", select = c("sleep_duration", "predicted_mexico"))

describe(mexico$predicted_mexico)
table(mexico$predicted_mexico)
describe(mexico$sleep_duration)

mexico_plot <- ggplot(mexico, aes(x=sleep_duration, y=predicted_mexico)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(1, 2.13, 4.1, 6.07, 8.04, 10.01, 11.98, 13.95, 15.92, 17.89, 20.17),
                     labels = c("1\n(Min.)", "2.13\n(μ-3SD)", "4.1\n(μ-2SD)", "6.07\n(μ-1SD)", "8.04\n(μ)", "10.01\n(μ+1SD)", "11.95\n(μ+2SD)", "13.97\n(μ+3SD)", "15.92\n(μ+4SD)", "17.89\n(μ+5SD)", "20.17\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.60, 0.9), breaks = seq(0.60, 0.9, by = .05)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
mexico_plot

mexico_plot <- ggExtra::ggMarginal(mexico_plot, type = "boxplot")

mexico_plot

issp_2007$predicted_new_zealand <- predict(poly.issp.m1, 
                                           newdata=transform(issp_2007, 
                                                             urban_rural_area_beta = 0,
                                                             age_beta  = 0, 
                                                             degree_beta  = 0,
                                                             income_beta   = 0,
                                                             sex_beta   = 0,
                                                             political_interest_beta   = 0,
                                                             religiosity_beta = 0,
                                                             country="New Zealand",
                                                             local_election=0,
                                                             mid_field_season=0,
                                                             day_off=0), type = "response", allow.new.levels=T)

new_zealand <- subset(issp_2007, subset = country == "New Zealand", select = c("sleep_duration", "predicted_new_zealand"))

describe(new_zealand$predicted_new_zealand)
table(new_zealand$predicted_new_zealand)
describe(new_zealand$sleep_duration)

new_zealand_plot <- ggplot(new_zealand, aes(x=sleep_duration, y=predicted_new_zealand)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(0.72, 3.03, 6.79, 8.67, 10.55, 14.31, 18.07, 22.5),
                     labels = c("0.72\n(Min.)", "3.03\n(μ-3SD)",  "6.79\n(μ-1SD)", "8.67\n(μ)", "10.55\n(μ+1SD)",  "14.31\n(μ+3SD)", "18.07\n(μ+5SD)" ,  "22.5\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.86, 0.94), breaks = seq(0.86, 0.94, by = 0.01)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
new_zealand_plot

new_zealand_plot <- ggExtra::ggMarginal(new_zealand_plot, type = "boxplot")

new_zealand_plot

issp_2007$predicted_philippines <- predict(poly.issp.m1, 
                                           newdata=transform(issp_2007, 
                                                             urban_rural_area_beta = 0,
                                                             age_beta  = 0, 
                                                             degree_beta  = 0,
                                                             income_beta   = 0,
                                                             sex_beta   = 0,
                                                             political_interest_beta   = 0,
                                                             religiosity_beta = 0,
                                                             country="the Philippines",
                                                             local_election=0,
                                                             mid_field_season=0,
                                                             day_off=0), type = "response", allow.new.levels=T)

philippines <- subset(issp_2007, subset = country == "the Philippines", select = c("sleep_duration", "predicted_philippines"))

describe(philippines$predicted_philippines)
table(philippines$predicted_philippines)
describe(philippines$sleep_duration)

philippines_plot <- ggplot(philippines, aes(x=sleep_duration, y=predicted_philippines)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(1, 3.6, 5.98, 8.36, 10.74, 13.12, 15.55, 17.88, 20.26 , 22),
                     labels = c("1\n(Min.)", "3.6\n(μ-2SD)", "5.98\n(μ-1SD)", "8.36\n(μ)", "10.74\n(μ+1SD)", "13.12\n(μ+2SD)",  "15.55\n(μ+3SD)",  "17.88\n(μ+4SD)",   "20.26\n(μ+5SD)",  "22\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.70, 0.92), breaks = seq(0.70, 0.92, by = .02)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
philippines_plot

philippines_plot <- ggExtra::ggMarginal(philippines_plot, type = "boxplot")

philippines_plot

issp_2007$predicted_south_korea <- predict(poly.issp.m1, 
                                           newdata=transform(issp_2007, 
                                                             urban_rural_area_beta = 0,
                                                             age_beta  = 0, 
                                                             degree_beta  = 0,
                                                             income_beta   = 0,
                                                             sex_beta   = 0,
                                                             political_interest_beta   = 0,
                                                             religiosity_beta = 0,
                                                             country="South Korea",
                                                             local_election=0,
                                                             mid_field_season=0,
                                                             day_off=0), type = "response", allow.new.levels=T)

south_korea <- subset(issp_2007, subset = country == "South Korea", select = c("sleep_duration", "predicted_south_korea"))

describe(south_korea$predicted_south_korea)
table(south_korea$predicted_south_korea)
describe(south_korea$sleep_duration)

south_korea_plot <- ggplot(south_korea, aes(x=sleep_duration, y=predicted_south_korea)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(0, 1.81, 4.54, 7.27, 10, 12.73, 15.46, 18.19, 22),
                     labels = c("0\n(Min.)", "1.81\n(μ-2SD)", "4.54\n(μ-1SD)", "7.27\n(μ)", "10\n(μ+1SD)", "12.73\n(μ+2SD)","15.46\n(μ+3SD)", "18.19\n(μ+4SD)",  "22\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.55, 0.90), breaks = seq(0.55, 0.90, by = 0.05)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
south_korea_plot

south_korea_plot <- ggExtra::ggMarginal(south_korea_plot, type = "boxplot")

south_korea_plot

issp_2007$predicted_netherlands <- predict(poly.issp.m1, 
                                           newdata=transform(issp_2007, 
                                                             urban_rural_area_beta = 0,
                                                             age_beta  = 0, 
                                                             degree_beta  = 0,
                                                             income_beta   = 0,
                                                             sex_beta   = 0,
                                                             political_interest_beta   = 0,
                                                             religiosity_beta = 0,
                                                             country="the Netherlands",
                                                             local_election=0,
                                                             mid_field_season=0,
                                                             day_off=0), type = "response", allow.new.levels=T)

netherlands <- subset(issp_2007, subset = country == "the Netherlands", select = c("sleep_duration", "predicted_netherlands"))

describe(netherlands$predicted_netherlands)
table(netherlands$predicted_netherlands)
describe(netherlands$sleep_duration)

netherlands_plot <- ggplot(netherlands, aes(x=sleep_duration, y=predicted_netherlands)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(1, 3.79, 6.19, 8.59, 10.99, 13.39, 15.79, 18.19, 22.5),
                     labels = c("1\n(Min.)", "3.79\n(μ-2SD)", "6.19\n(μ-1SD)", "8.59\n(μ)", "10.99\n(μ+1SD)", "13.39\n(μ+2SD)","15.79\n(μ+3SD)", "18.19\n(μ+4SD)",  "22.5\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.89, 0.92), breaks = seq(0.89, 0.92, by = 0.005)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
netherlands_plot

netherlands_plot <- ggExtra::ggMarginal(netherlands_plot, type = "boxplot")

netherlands_plot

issp_2007$predicted_russia <- predict(poly.issp.m1, 
                                      newdata=transform(issp_2007, 
                                                        urban_rural_area_beta = 0,
                                                        age_beta  = 0, 
                                                        degree_beta  = 0,
                                                        income_beta   = 0,
                                                        sex_beta   = 0,
                                                        political_interest_beta   = 0,
                                                        religiosity_beta = 0,
                                                        country="Russia",
                                                        local_election=0,
                                                        mid_field_season=0,
                                                        day_off=0), type = "response", allow.new.levels=T)

russia <- subset(issp_2007, subset = country == "Russia", select = c("sleep_duration", "predicted_russia"))

describe(russia$predicted_russia)
table(russia$predicted_russia)
describe(russia$sleep_duration)

russia_plot <- ggplot(russia, aes(x=sleep_duration, y=predicted_russia)) + 
  geom_point()  +
  stat_smooth(method='lm', formula = y ~ x, size = 0.5, level = 0.95, fullrange = TRUE, color = "black", linetype="longdash") +
  stat_smooth(method='lm', formula = y ~ poly(x,2), size = 0.5, level = 0.95, fullrange = TRUE, color = "black") +
  theme_classic() + ggpubr::grids(linetype = "dashed") +
  theme(axis.text.x=element_text(colour="black", size = 12), axis.text.y=element_text(colour="black", size = 12)) +
  theme(axis.title=element_text(size=14, colour="black", face="bold"), legend.text=element_text(size=14), legend.title=element_text(size=14)) +
  scale_x_continuous(breaks=c(0, 1.23, 5, 8.77, 12.54, 16.31, 20.08, 23.5),
                     labels = c("0\n(Min.)", "1.23\n(μ-2SD)", "5\n(μ-1SD)", "8.77\n(μ)", "12.54\n(μ+1SD)", "16.31\n(μ+2SD)", "20.08\n(μ+3SD)",  "23.5\n(Max.)")) + scale_y_continuous(breaks=number_ticks(10)) +
  scale_y_continuous(labels = function(x) paste0(x * 100, '%'), limits = c(0.95, 1), breaks = seq(0.95, 1, by = 0.01)) +
  labs(    x = "\n Sleep Duration (in Hours)\n", 
           y = "\n Fitted Values of Turnout \n",
           title = "")
russia_plot

russia_plot <- ggExtra::ggMarginal(russia_plot, type = "boxplot")

russia_plot

poly_figures <- ggpubr::ggarrange(finland_plot,  greece_plot, 
                                  ireland_plot, mexico_plot, netherlands_plot, new_zealand_plot,
                                  philippines_plot, russia_plot, south_korea_plot, 
                                  ncol=3, nrow = 3, labels = c("A: Finland", "B: Greece",
                                                               "C: Ireland", "D: Mexico", 
                                                               "E: the Netherlands",
                                                               "F: New Zealand",
                                                               "G: the Philippines", "H: Russia", "I: South Korea"))
poly_figures

ggsave(filename = "poly_figures_with_excluded_16.jpg",plot = poly_figures, width = 28, height = 16, dpi = 1000)
